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1.0 INTRODUCTION 

This document describes the Space Shuttle onboard navigation systems character- 
istics. Its purpose is to serve as a standard source of equations and numerical 
data for use in error analyses and mission simulations related to Space Shuttle 
development . 

The sensor characteristics described in this document are used for Shuttle 
onboard navigation performance assessment. Complete models are not always used 
for such studies, depending on the analyses to be performed, the capabilities of 
the computer programs, and the availability of computer resources. 

This document is not a control document . Controlling documents include- the Shut- 
tle Operational Data Book (refs. 1-1 and 1-2), the Coordinate Systems for the 
Space Shuttle Program (ref, 1-3), and the Space Shuttle Astrodynamical Constants 
(ref. 1-4) software configuration controlled input loads and trajectory design 
input data baseline . 

This document will be revised periodically. Suggestions for these revisions 
should be brought to the attention of R. T. Savely, FM8, at NASA/ JSC. 
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2 . 0 INE RTIAL MEASUREMENT UNIT 

2 . 1 GENERAL DESCRIPTION 

The inertial measurement unit (IMU) used on the Space Shuttle provides an 
inertially-f ixed stable platform. Mounted on the platfora are accelerometers 
that measure nongravitational accelerations and rate-integrating gyroscopes, 
which drive the gimbals. The platform is suspended by four gimbals to provide 
^TT-attitude rotational isolation of the platform from vehicle motion. The 
gimbals are read by girabal angle resolvers, which yield body attitude informa- 
tion with respect to inertial space (fig. 2-1). The outputs of the IMU 
(inertial accelerations and body attitude with respect to inertial space) are 
used by navigation and guidance for estimating position and velocity. 

The alineraent of the platform to a desired orientation is performed by first 
determining the present orientation and then accurately torquing the platform to 
the desired alinement. In-flight alinement is accomplished by optical sightings 
on stars using either the star tracker (sec. 5.0) or crew optical alinement 
sight (COAS) (sec. 6.0). Prelaunch alineraent is performed by gyro compassing/ 
leveling in which the gyros and accelerometers are used to orient the platform 
with respect to Earth rate and local gravity vectors. 
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VERTICAL GYRO 

Yg=Y Gyro I, A. 

Zq=Z Gyro I, A. 

AZIMUTH GYRO 

RG=Redundant Gyro I. A. 

Xj^=X Accel Sens . Axis 

DUAL AXIS ACCEL. 

^A“^ Accel Sens . Axis 

Zjj^=Z Accel Sens. Axis SINGLE AXIS ACCEL. 

T Gimbal Torquer 

FX=Gimbal Angle Resolver (1x/8x) 

CR=Gyro Coordinate Resolver 
RX1=Single Speed (8x) Resolver (Redundant Roll) 


Figure 2-1.- Inertial measurement unit component diagram. 
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■ ORIGINAL mm 

OF POOR QUALITY 79FH5 

2.2 SIMULATION APPROACH 

The IMU provides a measurement of the nongravitational acceleration applied to 
the IMU in the form of delta-velocity outputs, AVp, from a set of three 
orthogonally-mounted accelerometers. A measurement of the vehicle attitude is 
provided by the platform gimbal angles, 0, which define the orientation of the 
vehicle with respect to the inertial orientation of the platform'. 


An environment simulation provides the true vehicle acceleration, in a 

reference inertial coordinate system and a set of Euler angles, which de- 
fine the vehicle attitude relative to the reference system. Vehicle rate data 
may also be supplied by the environment, but we have chosen to derive the body 
rate within the IMU model. 


Figure 2-2 presents a block diagram for an IMU simulation. The functions of 
each block are as follows: 


a. Block 1 - computes the vehicle orientation matrix Tgj. 

Tbi = F (i) 


where F (^) is the transformation matrix that represents the rotation from 
reference inertial space to body coordinates defined by the Euler angle set ] 


b. Block 2 - the inertial acceleration, aj, is transformed to body 
coordinates . 


^ = '^BlSiI 

where Tgj is computed in block 1 . 
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Figure 2-2,- Inertial measurement unit simulation functional block diagram. 















c. Block 3 “ computes the average body attitude rate vector and the average 

body attitude acceleration vector over the time Interval At. The change in 
body attitude is given by 

R = (TBiJt - At 

The average body attitude rate vector, is computed for the eigenaxls ro- 

tation defined by R , 



assuming small attitude changes, sin ( [w|At)=^ jw|At. 

The average attitude acceleration vector is given by 

d. Block 4 - computes the sensed' acceleration, 3g, by adding rotational terms. 

-s “ ^B ^’^B ^ ^-B ^ + “b ^ - 

where ag was calculated in block 2, Wg and Kg were calculated in 
block 3, and r is the position vector of the IMU with respect to the center 
of gravity of the vehicle in body coordinates. 

e. Block 5 - expresses the sensed acceleration in ideal platform coordinates. 

ap = (tpi)(tbi)'^ ag 
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where ag was calculated in block 4, Tgj was calculated in block 1, and 
Tpj is the- transformation from reference inertial coordinates to ideal plat- 
form coordinates. 

f. Block 6 - computes the platfom drift vector, using the gyro error model 
equations presented in section 2.3.1. 

g. Block 7 - computes the platform initial misalinement error vector-, for 

either in-flight or prelaunch alinement using the equations presented in sec- 
tion 2.3*2. 

h. Block 8 - updates the platform misalinement matrix, M. The platform 
misalinement matrix is the rotation from ideal platform to actual platform 
coordinates. The change in M due to gyro drift, over time interval At 
is given by 

«t = R^t - At 

where ,R is the transformation matrix for the eigenaxis rotation about the 
drift vector d, assuming d is constant over the time interval At. 

R = I + (1-cos 6)u2/dm2 - (sin5)U/dm. 
where (1^ is the drift rate magnitude, | d| , 

I is the identity matrix 

S = dju At (eigenaxis incremental drift angle) 
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The initial misalinement matrix is formed from the misalinemerit error 
vector, computed in block 7 as follows: 

Mq = I + ( 1-gos4*)V^/4>^’- (sin<l5)V/<J> 


where (j) is the misalinement vector magnitude, [^q| 




' 0 

"'i'oz 

- 6 - 

O 

V 

= 

1 'I’oz 

0 

■<i>ox 



>> 

0 

1 

^ox 

0 


Block 9 - transforms the sensed acceleration from reference platform coordi- 
nates to actual platform coordinates. 



where a„ was computed in block 5 and M was computed in block 8. 

— P 

Block 10 - generates the simulated sensed delta- velocity vector AVp . 
AV p = (a 'p + At 

Where AVp is the accelerometer output before quantization 
Aa is the accelerometer error (sec. 2.3-3) 

At is the time interval between accelerometer samples 
^p is computed by block 9 and -assumed to be constant over time 
interval At. 
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Nondestructive quantization of the accelerometer data is performed in the 
following manner: 

It = it-it + «p 

= q(It) - q (Vp.ip) 

where is the accumulated delta-velocity vector at time t, generated by 

summing the nonquantized accelerometer outputs over the simulatio! 
interval 

q() is the quantization operator that truncates the argument to the 
least significant bit value for the specific accelerometer ■ being 
simulated 

AV p is the quantized accelerometer output vector. 

. Block 11 - computes vehicle orientation relative to the actual platform. 

TbP = Tbi TpiV 

where M is the misalinement matrix computed in block 8 
Tgj is the rotation matrix computed in block 1 

Tpj is the transformation from reference inertial to ideal platform 
coordinates . 

1. Block 12 - computes a set of gimbal angles, 0, from the IMU case orienta- 
tion vdth respect to the actual platform: 

■^CP = '‘’CN "^NB "^BP 
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where Tgp was calculated in block 11, Tjjg is the transformation from the 
body to navigation base coordinates (10.6° nav. base tilt): 


0 +sin 10.6° 

1 0 

0 cos 10.6° , 

J f 

and Tqjj is the IMU case mounting misalinement with respect to the navigation 
base matrix given by the equations in section 2.3.4. 


'NB - 


cos 10.6° 
0 

-sin 10.6° 


The platform is stabilized in an inertial orientation by the gyro output 
signals that drive torque motors on the azimuth, inner roll, and pitch 
gimbals. The outer roll gimbal receives its control signal from the 
inner roll resolver (the pick-off mounted between the inner roll and 
pitch gimbals). This control loop maintains perpendicularity between 
the inner roll and pitch gimbals, and the inner roll gimbal angle at 
zero degrees. 

The gimbal angle sequence defining the rotation from platform coordinates 
to vehicle coordinates is azimuth (®^), inner roll (5), pitch (®p), and 
outer roll (®jj). Since the outer roll gimbal control loop maintains • 
the inner roll gimbal angle (<5) at zero degrees, the vehicle attitude 
relative to the platform is defined by the Euler angle sequence 
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Then the case orientation matrix, T^p, is related to the gimbal angle set 
® by the equations given in section 2.3-5 that model the effects of 
gimbal nonorthogonalities . 


The gimbal angles, ©, are read by gimbal angle resolvers, and the result 
is quantized: 


0 = q (0 + ^) 


where 0 are the gimbal angles resulting from nonorthogonalities 
(sec. 2.3.5) 


A0 


are the gimbal angle resolver errors (sec. 2.3.6) 


q() is the quantization operator that truncates the argument to 
the least significant bit value for the specific resolver 
being simulated. 
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2.3 ERROR MODELS 

The mathematical models for accelerometers and gyros presented below are single- 
degree-of- freedom models. Two-degrees- of-freedom models can be constructed from 
these single-axis models by appropriate choice of error parameter values to 

model any cross -coupling or correlation exhibited by the two-degrees-of- freedom 
sensors . 

2.3.1 Gyro Error Model - Block 6 

Error in the gyro output can be modeled by bias drift, acceleration-sensitive 
(mass imbalance) drift, and acceleration-squared-sensitive (anisoelastic) drift. 
A mathematical error model for a rate- Integrating gyro is given by the equation: 

d = r’+^+h'^ag+ _ag'^ G ^g. 

where d is the total gyro drift rate 
r is the bias drift 
is the random drift 

is the acceleration vector in sensor input (i), spin 
(s), output (o) axis coordinates. Up is calculated in 
block 5. H is the transformation matrix from reference 
platform coordinates to sensor coordinates for the particu- 
lar gyro. 

is a vector of g-sensltive drift coefficients for accelera- 
tion along the input, spin, and output axes. 
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G = 


0 


0 


Sis 

®io 

®ss 

®so 

0 

Sqo 


is a matrix of g-squared sensitive drift coefficients 
for products of accelerations along the input, spin, 
and output axes. 


2.3.2 Platform Initial Misalinement - Block 7 

The generation of is dependent on the alinement technique and the sensors 

used in the alinement measurements. The simulation of an in-flight alinement 
using the star tracker (or COAS) and a prelaunch alinement using gyrocompassing 
and accelerometer leveling are presented here. 

2 . 3 . 2 . 1 In- flight Alinement 

The outputs from the star tracker Csec. 5 . 0 ) and the COAS (sec. 6 . 0 ) are 
measured line-of-sight (LOS) vectors to a pair of selected stars in actual plat- 
form coordinates. These LOS vectors are combined with the stored LOS vectors in 
reference inertial coordinates to form the transformation matrix from actual 
platform to reference inertial coordinates as follows (ref. 2 - 1 ). 

Let I31 and I32 be the stored LOS unit vectors in reference inertial 
coordinates. Define 


= Isi 

12 = UNIT (13^ X I32) 

13 = UNIT (Ii X I2) 

where UNIT() is an operator that normalizes a vector. The transformation from 
star coordinates to reference inertial coordinates is given by 

’■is = <Ii = Ja = I3> 
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Let Igi and Ig2 Le the measured LOS unit vectors in actual platform 
coordinates. Define 

h - ^si 

I2 = UNIT (Ig^ X Ig2) 

= UNIT (r, X I2) 

The transformation from actual platform to star coordinates is given by: 

"^SP “ • -2 * ^3^ 

Thus, the transformation from actual platform to reference inertial coordinates 
is 

"^IP = ^IS ”^3? ’ 

and the misalinement matrix, is given by 

^o = Tjp 

where Tjp is the transformation from ideal platform to reference inertial coor- 
dinates used in block 5. is computed from the off-diagonal elements of 

as follows. 


% (2,3) 

io = I -% (f.3) 


% ( 1 , 2 ) 


If the misalinements are too large the crew has several options for realinement 
and calibration of the IMU's. The nominal platform misalinement errors follow- 
ing in-flight star tracker realinement are presented in section 2.4.3. 
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2. 3. 2. 2 Prelaunch Alinement 

For prelaunch alinement the IMU’s are driven to a specified orientation, 
the desired platform orientation with respect to the navigation base. The 
specified orientation for IMU-1 is chosen such that IMU-1 X-platform axis is 
along the Up direction and tilted 0.4 degrees toward South, Y is pointed East, 
and 2 is pointed North and tilted 0.4 degrees above the local horizontal plane 
when the Shuttle is at the launch tower. IMU-2 and IMU-3 are related to IMU-1 
by the skewing matrix, Mpjp.], which describes the transformation from IMU-1 to 
IMU-2 and IMU-1 to IMU-3 for J = 2 and J = 3, respectively. P^jpi for launch 
is given in section 2.4.3. Also listed are the gimbal angles that the three 
IMU's will nominally display at platform release. 

Use of gyrocompassing and accelerometer leveling determine the actual platform 
orientation with respect to the North-West-Up coordinates at the time of platform 
release, %wuP* error inherent in this transformation is due to the 

accelerometers and gyros that are used, to determine the local-vertical and the 
Earth- rate vectors. Knowledge of the time of platform release, the Greenwich 
hour angle, and the geodetic latitude and longitude will yield the transformation 
from North-West-Up coordinates at the launch site to the reference inertial 
coordinates, The transformation from actual platform to reference iner- 

tial can now be calculated as 

"^IP = '^’iNWU %WUP 

Comparison of this transformation with the ideal platform to reference inertial 
transformation yields the misalinement matrix, M^, 

^o = ^IP^ Tjp = Th^^p'^ TfiwuP 
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Since accelerometer and gyro outputs are used to define the platform orientation 
with respect to the North-West-Up coordinates, the initial alinement errors, 
which are used to define , are highly correlated to the instrument errors . 


where 


<j) „ fd) (j) <j) 1 

_o - '■^on ow ou'' 

*^ou = ^u + ^U^GCA 

‘f^on = ‘i’N + 

“^ow = ‘^W 


T 


^u’ ^N» zero-mean random errors with a Gaussian distribution. Section 

2 . 4.3 gives the 1*^ values, dy is defined below. tQQ^ is the time from comple- 
tion of gyro-compass alinement to platform release (nominally 10 minutes). Zjj 
and are the north and west components of the accelerometer instrument 

biases (sec. 2.3.3) divided by the magnitude of the local gravity vector, g. 



From pJ.atform release to 12 seconds before lift-off, fine-tune gyro torquing, 
which' compensates for the known acceleration environment, is performed. The 
drift for this time period affects IMU misalinements as 

io = io d^t 

where t is the time from platform release and 

^ is the fine-tune gyro torquing compensation drift rate in North- 
West-Up coordinates 
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do = (dji dw du)*^ 

where icr values for djj and dy are given in section 2.4.3. The drift rate 
about the azimuth does not have the fine-tune compensation and is therefore 
given by the projection of the platform drift rate vector, calculated in sec- 
tion 2 . 3 . If into the azimuth (up) direction via Tj^^p. After lift-off minus 12 
seconds the platform misalinement is updated as shown in block 8. 


2 . 3.3 Accelerometer Error Model - Block 10 

Errors in the accelerometer output can be modeled by bias shift, scale factor 
error, input axis misalinement, and scale factor nonlinearity. A mathematical 
model for the acceleration measurement error is given by the equation 

Aa = b + h + k^^ + m (ai)^ + kalail 


where Aa is the acceleration measurement error 
b is the accelerometer bias 
Ti is the accelerometer noise 


k 



is a vector of scale factor error (k^) and input axis 
misalinements (Yj, Yk)* 




is the acceleration vector in sensor (i,j,k) axis 
coordinates. 3p is calculated in block 9- H is 
the transformation matrix from actual platform to 
sensor coordinates for the particular accelerometer. 
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is the acceleration component along- the sensor input axis, 

are the acceleration components along orthogonal axes perpendicu- 
lar to the input axis, 
m is the scale factor nonlinearity. 

is the scale factor asymmetry error. 

The error model parameters are assumed to be independent random variables with 
zero-mean Gaussian distributions. 

2.3.^ IMP Case to Navigation Base Mounting Misalinement - Block 12 
The IMU case to navigation base mounting misalinement is a function of not only 
the accuracy of the case installation, but also the stiffness of the navigation 
base itself. Let C be a vector of the IMU case misalinements with respect to 
the navigation base. The transformation from navigation base to the IMU case is 
given by 


Tcu = I + (1 - cosC) U^/C^ •- (sinC) U /C 

where I is the 3x3 identity matrix 

C = |c| , the magnitude of the IMU case misalinement vector 



-Cy 0 


The error model parameters are assumed to be independent random variables with 
•zero-mean Gaussian distributions. 
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2.3-5 Gimbal Angle Nonorthogonality Errors.- Block 12 

The following discussion is derived entirely from reference 2-2. Let 0^, 0p, 
0p represent the inner, middle, and outer gimbal angles; 6, the inner roll 
gimbal angle 

Y = the small angle misalinement of the azimuth axis towards the inne: 
roll axis 

y = the small angle misalinement of the inner roll axis towards the 
pitch axis 

e = the small angle misalinement of the pitch axis towards the outer 
roll axis 

The outer roll axis of rotation is defined along the X ease axis. Any 
misalinement of the outer roll axis to the IMU case can be included in the IMU 
ease to navigation base misalinement (sec. 2.3-5). 

The transformation fr<xn actual platform to IMU case coordinates can be written 
as a series of seven single axis transformations 


where : 


R2(0) 


C0 

S0 

0 

-S0 

C0 

0 

0 

0 

1 
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and S0 = sin 0 and C0 = cos 0. 

Equating the above matrix, after multiplication of the seven rotation matrices 
using small angle approximations for Y, y, e, and 6 and neglecting sec- 
ond-order terms, with the matrix 

^CP = '^BP = N 

calculated in block 12, leads to the following iterative solution. 

A 

= tan"'’ ( N(1,2)-£C0A-^C®AS®p3/tN(1 ,1)-e SO^ - Ss0AS0p] - U 

/S 

0p = sin-’’ (-N(1,3))-Y 

0R = tan"'’[N(2,3)-es0pC0R-5c0R]/(N(3,3)+es©pS0R+6sQR) 

As an initial estimate of the real girabal angles, use 

0A = tan"’’ N(1,2)/N(1 ,1) 

0p = sin""' -N(1,3) - Y 
0R = tan""' N(2,3)/N(3,3) 

6=0 
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The error model parameters are assumed to be independent random variables with 
nonzero-mean Gaussian distributions. The azimuth and pitch nonorthogonalities 
(Y and p) are compensated for in the IMG SOP by adding Y to the pitch 
resolver readout (0p) and by adding y to the yaw resolver readout (0^.). The 
outer roll nonorthogonality (g) is measured in Hanger Cal B and is compensated 
for in the attitude processor onboard software module. Values for the range of 
the means and the errors in calibration are given in section 2.4.6. 

2 . 3.6 Resolver Error Model - Block 12 

Gimbal angle resolver errors can be categorized as those caused by bias shift 
and sinusoidal resolver errors. A mathematical model for the gimbal angle 
resolver error is given by the equation: . 


A0 = Bg + Tig +^2^Aq^ sin (n0 + <|)jj) 


where A0 

B0 

^0 

A0n 

4>n 

A 

0 


is the total resolver error 

is the resolver bias 

is the resolver random noise 

is the sinusoidal bias for multiplicative speed n 
is the random phase error for multiplicative speed n 

A A A 

is the gimbal angle component (0^, 6, 0p, or 0p) resulting 

from nonorthogonalities, as computed in section 2.3.5. 


The error model parameters are assumed to be random independent variables with 
zero-mean Gaussian distributions, except for the phase angle errors, vdiich 

have zero-mean uniform distributions. 
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2.4 ERROR SOURCE VALUES 


The Space Shuttle contains three Singer-Kearfott IMU's. Each IMU contains two 
two-degrees -of- freedom Gyroflex gyros and one single-axis accelerometer and one 
dual-axis accelerometer. 


2.4.1 IMU Locations 

The location of the three IMU's are given below in Orbiter structural body coor- 
dinates with units of inches 


IMU 



Zo 

1 

419.187 

-13.75 

422.0 

2 

419.187 

0 

422.0 

3 

419.187 

13.75 

422.0 


To determine the locations in body coordinates, the center of mass of the 
Orbiter for the particular flight must be expressed in Orbiter structural body 
coordinates and subtracted from the IMU locations listed above. The following 
transformation matrix, ^BO' applied to the resultant difference vector will 
give the IMU locations in body coordinates. 


-1 0 

0 1 

0 0 


0.015 deg/hr for 17 hours after 
preflight alinement and calibra- 


2.4.2 Gyro Error Values (in) 
Bias drift (r) - all axes 


0 

0 

-1 
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g- sensitive drift - input and output axes 


(hi,h„) 


- spin axis (hg) 


0.022 deg/hr for 17 hours after 
each onorbit calibration 
0.035 deg/hr thereafter 
0.025 deg/hr/g 

0 deg/hr/g 


g^-sensitive drift - off-diagonal terms 

Sio’ ®so^ 

- diagonal terms 
Sqo’ ®ss^ 


0.025 deg/hr/g^ 
2 

0 deg/hr/g 


Transformation matrix for X-axis gyro 


% = 


Transformation matrix for Y-axis gyro 


1 0 0 

0 0-1 

0 -1 . 0 


% = 


0 1 0 

1 0 0 

0 0-1 


Transformation matrix for Z-axis gyro 


Hz = 


0 

-sin 12° 
cos 12° 


0 

cos 12° 
sin 12° 


1 

0 

0 
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2.4.3 Initial Mlsalinement Error Values (la) 


2.4.3.-I In-flight Alinement in Navigation Base Coordinates 


ox 

71.2 

sec 


72.4 


oy 

sec 


71.1 


oz 

sec 


The following is the RSS of: 

Star tracker bias pointing error, all axes 


Star tracker random error (21 measurement average), 
all axes 

Star tracker to IMU case misalinement , X-, Y-axis 

Z-axis 

IMU case to platform bias error X-axis 


Y-axis 

Z-axis 


60/ sec 
10.6/y'5T sec 


32.2 

29.2 
47.1 
49.0 
49.0 


sec 

sec 

sec 

sec 

sec 


The following is the RSS of: 

Outer roll to case gimbal nonorthogonality, Y- and Z-axes 
Resolver bias, Y- and Z-axes 
X-axis 

Outer roll to pitch gimbal nonorthogonality, Y- and Z-axes 


20 sec 

49.0 sec 

42. 1 s^ 
3O//2 s^ 


The star tracker to IMU case misalinement is the RSS of: 
Star tracker to mounting pads, all axes 
Mounting pad to nav. base, all axes 
Thermal nav. base bending, X- and Y-axes. 

Nav. base to mounting pads, all axes 
Mounting pad to IMU case, all axes 


20 sec 
5 sec 
41/3 sec 


5 sec 
20 sec 
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4>U 

fine-tune gyro torquing compensation drift rates 
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8.7 sec 
60 sec 

.0094 sec/sec 


Completion of gyrocorapassing nominally occurs 600 seconds before platform 
release. Accelerometer leveling is completed at platform release, and lift-off 
nominally occurs 20 minutes after platform release (not including the possibil- 
ity of a 10-minute and/or a 20-minute hold) . 


IMU gimbal angles at platform release 


Roll ( deg ) 
0.0000 
334.4069 
296.2381 


Pitch (deg) 

349.0000 

44.3269 

25.5850 


Yaw (deg) 
0.0000 
309.7792 
350.4931 


(IMU1) = 


(IM02) = 


0.9816272 

0 

-0.1908089 

0 

1 

0 

.1908089 

0 

.9816272 

0.645-1763 

-0.2079225 

0.7351977 

-.3090204 

.8090255 

.4999842 

-.6987517 

-.5497691 

.4577118 
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-.8090132 

-.4318494 


.5000085 

-.1489717 


0.3364485 
-.3090127 
.8895581 , 
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TP2P1 = Tdn(IMU2) Tdn(IMU1)^;Tp3pi = Tdn(IM03) Tdn(IMUI)' 


2.4.4 Accelerometer Error Values (icr) 

Bias (b) - all axes 

Noise (n) - all axes 

Scale factor (k^^) - all axes 


Input axis alinement (Yj,Yj^) - all axes 
Scale factor nonlinearity 
Scale factor asymmetry 
Quantization 


Transformation matrix for X-axis accelerometer 


H 


X 


1 

0 

0 


0 0 
1 0 
0 1 


Transformation matrix for 


Y-axis acceTerometer 


H 


y 



1 

0 

0 


0 

0 

1 


50 yg 
5 yg 

40 PPM for 17 hours 
after preflight aline- 
ment and calibration 
100 PPM thereafter 
15 

0 g"^ 

40 PPM 

1 cm/ sec 
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Transformation matrix for Z-axis accelerometer 


Hz 


0 

0 

1 


0 

1 

0 


1 

0 

0 . 


2.U.5 IMU Case Error Values (la) 

IMU case to navigation base misalinement (Cjj, Cy, Cy) 


These are the RSS of the following error sources 

I 

[ IMU case to mounting pads - all axes 20 

Mounting pads to nav. base - all axes 5 sec 


Nav. base bending (static loads) - X, Y axes 
Ascent 
Onorbit 
Entry 


5//2 sec 
0 sec 
4//2 s^ 


Nav. base bending (thermal) - X, Y axes 

Onorbit only ' 41/3/2 

Outer-roll to case gimbal nonorthogonality - 20 s^ 

Y, Z axes 


2.4.6 Gimbal Nonorthogonalities (y» u» e) 
Range of means (y, p, e) 

Error in calibration (e) - 1o 
Error in calibnation (y, p). - Ict 


+413 

,30 

0 


sec. 


sec 


sec 
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2.4,7 Resolver and A/D Converter Error Values (1g) 
Resolver bias (B 0 ) 

Random noise (t ]0 ) 

Sinousoidal bias 

First harmonic (Aq,| ) 

Eighth harmonic (A 0 g) 

Ninth harmonic (Agg) 

Sixteenth harmonic (A 0 .|g) 


30 see 
12 sec 


7.6 sec 
19.0 sec 
4.2 s^ 


20.0 sec 


Random phase error (4>n) -180°' < (|)jj < 180° 

Quantization 20.0 sec 




Of POOR 
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3.0 BODY MOUNTED SENSORS 

3.1 GENERAL DESCRIPTION 
Body mounted sensors consist of rate gyro assemblies (RGA) and accelerometer 
assemblies (AA). The RGA contains three rate gyros that are mounted to the ve- 
hicle frame with the input axes forming an orthogonal triad (fig. 3-1). The 
outputs from the RGA provide a measurement of vehicle attitude rate used to 
update the onboard estimate of vehicle inertial attitude. The AA contains nor- 
mal and lateral accelerometers mounted to the vehicle frame (fig. 3-1). The 
outputs from the AA provide a measurement of the nongravitational accelerations 
applied to the vehicle. The AA outputs are used by Flight Control to provide 
load relief during the high dynamic pressure region of the ascent phase and 
steering commands to the entry phase digital autopilot to conform to a selected 
g-profile trajectory. 

^C£D»^'PAQE NOT 



Figure 3-1.- Body mounted sensors orientation. 
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3.2 SIMJUTION APPROACH 

The body mounted sensors provide measurements of the normal and lateral 
nongravitational body accelerations, (n =1, 2), and measurements of the 
body attitude rates, (n = 1, 2, 3)- 

An environment simulation is assumed to provide- the true vehicle acceleration, 
aj, in a reference inertial coordinate system and a set of Euler angles, 3, 
which define the vehicle attitude relative to the reference system. Vehicle 
rate data may also be supplied by the environment, but the assumption made here 
is that the rate must be computed within the simulation. 

Figure 3.2 presents a block diagram for a body mounted sensors simulation. The 
functions of each block are as follows. 

a. Block 1 - computes the vehicle orientation matrix Tgj. 

Tgi = F (3) 

where F (3) is the transformation matrix that represents the rotation 
defined by the Euler angle set 3. 

b. Block 2 - the inertial reference acceleration, a^, is transformed to body 
coordinates . 

fg = ^ 

where Tgj is computed in block 1. 

c. Block 3 - computes the average body attitude rate vector and the average 
body attitude acceleration vector over the time interval At. The change in 


body attitude is given by 
R = (TbiJ\ -At 
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Figure 3-2.- Body mounted sensors simulation block diagram 
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The average body attitude rate vector, Ug, is computed for the eigenaxis ro- 
tation defined by R. 



assuming minor attitude changes, sin (|wgiAt) IwglAt. 

The average attitude acceleration vector is given by 
c?b = ([“b^t ■ ^“B^t -At^ 

d. Block 4 - computes noninertial acceleration by adding rotational terms. 

✓s. 

.5b “ .5b ^!^B ^ ^ + “B ^ r 

where ag was calculated in block 2, Wg and ctg were calculated in block 
3, and r Is the position vector of either the rate gyro assembly 
and the acceleration is denoted by or the accelerometer assembly 

(r^^^ and the acceleration is denoted by with respect to the. cen- 

ter of gravity of the vehicle. 

e. Block 5 - generates the simulated nongravitational acceleration sensed by 
the accelerometer assembly. 

a^i = q + Aa„) 

■aj = q (a^ + Aa^) 

where and are the respective z- and y-components 

of the body acceleration, ag^'^^ computed in block 4, Aajj and 
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Aa;^ the normal and lateral accelerometer errors (sec. 3.3-1) 'and 

q ( ) is the quantization operator that truncates the argument to the least 
significant bit value for the specific accelerometer being simulated. 

Block 6 - computes the simulated rate gyrO' output 

idB = ‘I + Awb) . 

where ^ are the body attitude rates computed in block 3 , Aug are 
the rate gyro errors (sec. 3.3.2), and q ( ) is the quantization operator 
similar to the one in block 5. 
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3.3 ERROR MODELS 
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3-3.1 Aocelerometer Assembly Error Model - Block 5 

Single-axis accelerometer errors can be categorized as those caused by bias, ra 
dom noise, scale factor error, and input axis misalinement. A mathematical 
model for the normal accelerometer is given by the equation 


Aar 






where Aa^j is the normal accelerometer measurement error 
bjj is the normal accelerometer bias 
is the normal accelerometer noise 


is a vector of longitudinal and lateral 
misalinements ("Ynx' "^ny^ scale factor 
error (k^.). 

'' (AA) ■ 

a-Q^ ■' zs the acceleration vector of the accelerometer assembly in body c 
ordinates (x, y, z - longitudinal, lateral, normal). 



The lateral accelerometer model is given by the equation 


Aai = 


where Aa^^ 


is the lateral accelerometer measurement error 

is. the lateral accelerometer bias 

is the lateral accelerometer noise 



is a vector of longitudinal and normal misalinements 
(Yix> YIz^ scale factor error (ky) 
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is the acceleration vector of the accelerometer assembly in body 
coordinates . 

The error model parameters are assumed to be independent -random variables with 
zero-mean Gaussian distribution. 


3 . 3.2 Rate Gyro Assembly Error Model - Block 6 

Single-degree-of- freedom rate gyro errors can be categorized as those caused by 
bias drift, random drift, linear acceleration sensitive drift, angular accelera- 
tion sensitive drift, scale factor error, and input axis misalinements . A mathe- 
matical model for the rate gyro measurement error is given by the equation 


AWb = 

where Awg 


r 


n 

h 

aB'”") 


g 

“b 


r . n . + g ag -f k'^tiOg 

is the total rate measurement error for a given axis (x,y,z - roll, 
pitch, yaw) 

is the bias drift rate 

is the random noise drift rate 

is the g-sensitive drift rate 

is the linear acceleration component for the rate gyro assembly 
along the given axis 

is the angular acceleration sensitive drift rate 

is the angular acceleration component along the given axis 


k 



is a vector of scale factor error (kj^) and input 
axis misalinements (Yj, Y,,) 
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^ is the body attitude rate vector 

The error model parameters are assumed to be independent random variables with 
zero-mean Gaussian distribution. 
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3.4 ERROR SOURCE VALUES 
3.4.1 Accelerometer Assembly 

The accelerometer assembly is composed of two Honeywell single-axis acceler- 
ometers (model No. 326) mounted orthogonal to each other. Four of these packages 
are then attached to the vehicle at the locations specified in section 3.4.3. 
References 3-1 through 3-3 provide the following error source values. 


Accelerometer errors (Ih) 



Bias - normal (bj^) 

8333 

yg 

lateral (bj^) 

5000 

yg 

Random noise - both axes rii) 

3333 

yg 

Scale factor - normal (k^) 

0.833 

percent 

- lateral (ky) 

1.0 percent 

Input axis misalinements (all axes) 
This is the resultant of: 

0.22° 


Input axis to case misalinement 

0.12° 

Case-to-body misalinement 


0.18° 

Quantization - normal 

7800 

yg 

- lateral 

2000 

yg 
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3.^.2 Rate Gyro Assembly 

The rate gyro assembly is composed of three Northrop single degree-of-freedom 
rate gyros (model No. G-6). Four assemblies are attached to the vehicle at the 
locations specified in section 3 •^•3* References 3-1 through 3-3 provide the 
following error source values. 

Rate Gyro Errors (la) 

Bias - all axes (r) .05 deg/sec 

Noise - yaw and pitch ('Tiy} Hg) .01667 deg/sec 

- roll (Ojj) .03333 deg/sec 

Linear acceleration sensitivity •- all axes (h) 0.01667 deg/sec/g 

' P 

Angular acceleration sensitivity - all axes (g) 0.001 deg/sec/deg/sec 

Scale factor - all axes (k) 1.7 percent 

Input axis misalinement ~ all axes 0.248° 

Quantization - pitch and yaw 0.039 deg/sec 

- roll 0.078 deg/sec 

3.4.3’ Body Mounted Sensor Locations 

The locations of the four accelerometer assemblies and the four rate gyro 
assemblies are given below in Orbiter structural body coordinates with units of 
inches . 
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Xo 

Xo 

Zo 

AA1 

382.8 

-6.8 

398.2 

AA2 

381 .7 

5.6 

406.4 

AA3 

381.7 

11.6 

406.4 

AA4 

382.8 

-12.8 

398.2 

RGA1 

1307 

-100 

353 

RGA2 

1307 

100 

353 

RGA3 

1307 

-13 

295 

RGA4 

1307 

-21 

101 


To determine the locations in body coordinates, the center of mass of the 

Orbiter for the particular flight must be expressed in Orbiter structural 

body coordinates and subtracted from the sensor locations listed above. 

The following transformation matrix, T, applied to the resultant difference 
! 

vector will give the sensor locations in body coordinates. 


T 


-1 

0 

0 


0 -0 
1 0 
0 -1 
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4.0 BAROMETRIC ALTIMETER 

4.1 GENERAL DESCRIPTION 

A barometric altimeter is a device that is used to determine the ambient 
atmospheric pressure around a vehicle. Use of an atmosphere model then allows 
pressure to be related to vehicle altitude through a pressure/altitude 
algorithm. The device itself only measures pressure, and the transformation to 
altitude is performed within the onboard software. Barometric altimeters are 
generally restricted to altitudes of less than 100 000 feet. 

4.2 SIMULATION APPROACH 

The basic data type for simulation of a barometric altimeter is the geometric al- 
titude. The geometric altitude is the altitude (h) above the ellipsoidal Earth 
(fig. 4-1). Figure 4-2 presents a block diagram for the barometric altimeter 
simulation. A functional description of each block is as follows. 

a. Block 1 - calculates the geometric altitude (h) for input to the baro error 
model block as follows (ref. 4-1). 

h = (1-RE(1-e)yz2E7+ (1-e)^ R^xy)(Z^EF+n-e)^R^XY)/'^Z^F+ 

vfliere 

Rg = equatorial radius of the Earth 
e = Earth flattening 
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Figure 4-1.- Altitude above an ellipsoidal Earth. 
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Figure 4-2.- Barometric altimeter simulation block diagram. 
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Zgp = Z - component of vehicle position in an Earth-fixed coordinate system 
= X^eF + 

b. Block 2 - computes a corrupted geometric altitude by modifying the output of 
block 1 with the following error sources (ref. H-2) . 

H — h + 

where 

AHgg = Bias error 

AHpE = Random error in the air data transducer assembly (ADTA) 

c . Block 3 - compensates the corrupted altitude (H) from block 2 by adding a 
correction term for nonstandard atmosphere error (AHgp) or subtracts 

a barometric pressure update error (AH^ggp), depending on whether the 
corrupted altitude is above or below a specified altitude. The following 
equations represent this procedure. 


If H > 

^^'^SPEC 

Then: 

Hq = H + AHgp 

Else: 

He = H - ■^LSBP 


where 

ALTgpEc = specified a^ltitude 

H = corrupted geometric altitude calculated in block 2 

AHgF = Nonstandard atmosphere error 
AHl 2 bp - barometric pressure update error 
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k.3 ERROR MODEL 

The environment baro altimeter error model presented in reference k-2 identifies 
the following four error sources. 

a. Bias error - a bias error due to the imperfect measurement of pressure by 
the barometer (ADTA). 

b. Barometric pressure update error - Uncertainty in the landing site baromet- 
ric pressure, which is used to update the" vehicle’s barometric pressure 
shortly before landing. 

c. Scale factor error - an error due to the difference in average pressure at 
a particular location not agreeing with the average pressures shown in the 
1962 Standard Atmosphere (ref. 4-4). 

d . Random error - A Gaussian noise inherent in the ADTA . 

The one- sigma values for the baro altimeter error model are presented in table 
4-1. The bias error, uncertainty in the landing site barometric pressure, and 
scale factor error comprise the correlated error, while the noise error com- 
prises the uncorrelated error that is added to the barometer's measurements. 

The bias error plot in table 4-1 was constructed from a simulation of an OFT-1 
trajectory to EAFB (Edwards Air Force Base). The uncertainty in the landing site 
barometric pressure and scale factor error are modeled as ECRV's (exponential 
correlated random variable) with time constants of t = “> and t = 1000 sec, re- 
spectively. The uncorrelated error is also modeled as an ECRV with a time con- 
stant of T = 0, vdilch implies a true random error. 
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The specified altitude (ALTgpgQ^ for adjusting the compensation errors in block 3 
is given as 5000 feet. 
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TABLE 4-1.- BAROMETRIC ALTIMETER ERROR MODEL (1 -sigma) 

Correlated error 

( 1 ) Bias error 



altitude (kft) 

(2) Barometric pressure update error • ■ 

AHl 3 bp = 20 ft (t = «>) vehicle is updated with the landing site 

pressure shortly before landing 

(3) Scale factor error related to density reference 
AHgF = 3 - 5 % of altitude (t = 1000 sec) 

using the 1962 standard atmosphere model in 
the filter 

(4) Instrument Noise AHpg =2.3 exp (h/ 23000) ft (x = 0) 
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5.0 STAR TRACKER 

5 . 1 GENERAL DESCRIPTION 

The Shuttle star tracker is a strapped- down, wide field-of-view image-dissector, 
electro-optical searching and tracking device. The star tracker is used to ob- 
tain precise angular measurements on stars or Sun- illuminated targets of small 
angular image diameter. There are two star trackers mounted on an extension of 
the Shuttle navigation base (fig. 5-1) . The -Z star tracker points in the ap- 
proximate direction of the Shuttle -Z body axis, while the -Y star tracker 
centerline is approximately 10.5 degrees forward of the Orbiter -Y body axis. 

The star tracker instrument consists of both the actual star tracker and a light 
shade. The light shade enables the tracking of targets that are relatively 
close to the Sun, Moon, or the sunlit Earth's horizon. The star tracker has the 
capability to search for, acquire, and track the 153 brightest stars or sunlit 
targets of 3 to -7 magnitude at any location within its 10 x 10 degree square 
field of view. 

The two star trackers, under onboard computer software control, will track stars 
for the purpose of IMU platform realinement. Since the IMU inertial platform 
drifts from a true inertial attitude its- orientation must be periodically deter- 
mined in order to make corrections. Using star tracker measurement data for two 
stars, and the known star positions, the IMU inertial platform orientation with 
respect to the mean of 1950 (M50) coordinate system can be computed using 
deterministic methods. The accumulated IMU error due to drift can then be 
removed by torquing the -IMU platform back to its original position. The star 
trackers can also be used to track a sunlit rendezvous vehicle. In this case 
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the star tracker measured angular position data is incorporated into the 
rendezvous navigation filter for use in relative state estimation. 

5.1.1 Star Tracker Operating Characteristics 

The Shuttle star tracker contains the following three modes of operation: 

a. Automatic scan mode 

b. Offset scan mode 

c. Self test mode 

In the automatic scan mode (the default mode when the star tracker is first 
turned on), the star tracker scans the entire 10- x 10-degree field of view 
from top to bottom repeatedly. The search scheme (fig. 5-2) is a left-to- 
right, right- to-left, top-down raster scan. If a target encountered in the 
search has a magnitude greater than the currently commanded sensitivity thresh- 
old setting, the target will be acquired and the track scan will begin. The 
track scan is a cruciform pattern centered on the position of the target in the 
field of view (fig. 5-2). If the target is moving the track scan will follow 
it, remaining centered on the target. During the track scan, the measured hori- 
zontal and vertical components of the position of the target, and the measured 
target magnitude, are output by the tracker at a 25-Hz rate. The target will 
continue to be tracked until either it leaves the star tracker field of view or 
a break track command is sent to the tracker. .In either case the star tracker 
will return to the auto scan mode and search the remainder of the field of view. 

When the star tracker is commanded to the offset scan mode it will search a re- 
duced 1- X 1-degnee field of view. centered at a position in the full field of 
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view defined by a pair of commanded offset coordinates. The reduced field-of- 
view search is also a top-down raster scan similar to the auto scan search. 

When a target is acquired in the offset mode it will be tracked until either a 
break track is commanded or the target moves outside of the full 10- x 10-degree 
field of view. In either case the star tracker will return to the originally 
commanded offset coordinates and resume the reduced field-of-view scan. 

The star tracker can also be commanded to a self- test mode. When this mode is 

I 

commanded the star tracker searches for and tracks an LED source that is re- 
flected into the star tracker field of view near -the edge via a small aluminized 
spot on the star tracker protective window. The measured offset angles and mag- 
nitude of the simulated star are compared with prestored values to determine 
whether the star tracker operational accuracy conforms with requirements. 

5-1.2 Star Tracker Operating Constraints 

The following list briefly summarizes the various star tracker specifications 
and operating constraints. 

a. Field of view dimensions 
Auto scan - 10 x 10 degrees 
Offset scan -1x1 degrees 

b. Field of view search time 
Auto scan - 10 seconds 
Offset scan - 1 second 

c. Target Intensity limits 
Minimum intensity - 3-0 magnitude 
Maximum intensity - 7.0 magnitude 
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d. Target image diameter must be no greater than 8 arc minutes. 

e. The star tracker does not have the capability to track either star of a 
pair that have an absolute magnitude difference greater than one and an an- 
gular separation' less than one-half degree. 

f. The star tracker will acquire and maintain accurate tracking of stars for 
Orbiter angular rates of up to 0.2 deg/ sec and will maintain lockon with 
reduced accuracy at rates of up to 0.5 deg/sec. 

g. Bright source avoidance angles' 

For acquisition and accurate tracking of stars or targets the minimum sepa- 
ration angles between the star tracker boresight axis and the three possi- 
ble sources of intense light are as follows: 



(1) 

Center of Sun 

0 

o 

CO 


(2) 

Center of Moon 

8“ 


(3) 

Earth’s horizon 

20° 

h. 

Sensitivity threshold 

levels 


The star tracker will not acquire stars or targets having intensities that 
are less than the currently commanded sensitivity threshold level. The 
four crew commandable threshold levels correspond to star magnitudes of 
3.0, 2.^, 2.0, and 1.0. 

i. Offset, command angle quantization 

The reduced field- of-view offset command angles have a granularity of one- 
third degree. 
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j . Warm-up time 

Star tracker performance requirements are met within 15 minutes after power 
on . 

k. Accuracy requirements 

(1) Angular data 30 arc sec noise (icr) 

30 arc set bias (Ict) 

(2) Intensity error +0.6 magnitude 

(a detailed angular data error budget is presented in sec. 5.4) 

5.1.3 Star Tracker Hardware/Software System Operation 

The Shuttle onboard star tracker computer software accumulates star tracker angu- 
lar measurement data for each star tracked and converts it into an average 
measured unit line-of- sight vector expressed in IMU platform coordinates. Two 
star line-of-sight vectors with adequate separation are required by the IMU on- 
orbit alinement software for the purpose of realinement of each IMU platform. 

Nominal star tracker operation assumes the execution of a vehicle attitude rate 
by the crew. As each star tracker field of view sweeps across the celestial 
sphere, the software selects stars of opportvinity for tracking. Stars are 
selected from a navigation star catalog based on an estimated star tracker line 
of sight. The estimated star tracker line of sight is computed based upon a 
reasonably well known IMU inertial platform attitude with respect to the celes- 
tial sphere. Prior to star selection software checks are performed to determine 
whether either star tracker is occulted by the Earth, Sun, or Moon. If it is 
determined that a star is visible in either star tracker field of view, the star 
tracker is commanded by the software to search a reduced 1- x 1-degree field of 
view centered at the estimated position of the star in the star tracker full 
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field of view. If a star is acquired a sequence of 21 measured star tracker 
output angles are accumulated at a 6.25-Hz rate and then used to compute the av- 
erage measured horizontal and vertical offset angles of the star. Simultaneouslj 
a sequence of 21 IMU gimbal angles are accumulated and averaged. Twenty- one 
samples are accumulated because this number of samples is the smallest number 
that yields the minimum average time lag between sensor output and the software 
read. The time lag is a result of the asynchronous operation of the star tracker 
and the software. The measured star tracker and IMU data are then converted to 
a unit line-of-sight vector in IMU coordinates by the following equations. 


a. First, compute the navigation base to IMU cluster transformation matrix 
using the following average IMU gimbal angles: 

AZ = average azimuth gimbal angle 
IR = average inner roll gimbal angle 
P = average pitch gimbal angle 
OR = average outer roll gimbal angle 
SAZ = sin AZ 
CAZ = cos AZ 


SIR = sin IR 
CIR = cos IR 
SP = sin P 
CP = cos P 


SOR = sin OR 


COR = cos OR 


(tnbpc] = 


CAZ 

SAZ 

0 


-SAZ 

CAZ 

0 


0 

0 

1 


10 0 


" CP 0 SP‘ 

0 CIR -SIR 


0 10 

0 SIR CIR 


-SP 0 CP 

- 


— ■■ 
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1 -DP 0 

DP 1 0 

0 0 1 


10 0 
0 COR -SOR 

0 SOR COR 


[tnbrl) 


DP = the nonorthogonality between the outer roll and pitch axes 
(constant) . 

(tnbrl] = navigation base to IMU roll axis transformation (constant) 


b. Compute the observed star line-of-sight vector in IMU platform coordinates 
as follows: 


s = (tnbpc) [tnest]"^ 

V = average measured 
H = average measured 
(tNBST] = navigation 

c. Correct the observed line of sight for stellar aberration 



' 

/-sin g\ 1 

COR = , 

V + s (teclmso) 

1 cos G j 1 

V ° /j 


TECLM50 = ecliptic to mean of 1950 coordinate transformation 
V = vehicle velocity vector in M50 coordinates 
G = mean longitude of the Sun at the current time 
S = Earth's orbital speed 
c. = speed of light 
S = unit (S - (tcMSo)'^ ^) 

(tcM5o] = present IMU cluster to M50 transformation matrix 


/-tan V \ 

UNIT ^ tan H j 

vertical offset angle 

horizontal offset angles 

base to star tracker transformation 


5-9 



79FM5 


A measured unit star line-of-sight vector in IMU platform coordinates is 
computed for each operating IMU. The star tracker software, in addition to 
accumulating star tracker measurements and computing the equivalent line-of- 
sight data, also performs data reasonableness checks for each star and optimum 
data selection operations. When adequate measurement data have been gathered 
for two stars an IMU realinement can be performed by the IMU onorbit alineraent 
software . 

\ 

For the purpose of rendezvous vehicle tracking initial acquisition of the 
rendezvous target is accomplished by a crew maneuver of the Orbiter to center 
the target in the selected star tracker field of view and execution of the star 
tracker software to command an offset mode search for the target. The onboard 
rendezvous navigation software then accepts average measured star tracker output 
angles from the star tracker software every 15 seconds for incorporation into 
the rendezvous navigation filter. 
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5.2 SIMULATION APPROACH 
The diagram in figure 5-3 illustrates the method used in the simulation of the 
star tracker hardware. It is assumed that the following eight inputs are avail- 
able frcm an environment simulation. 

a . [ TIB) = mean of 1950 to vehicle body transformation 

b. (TBN) = vehicle body to navigation base transformation 

c. (tNS) = ideal navigation base to star tracker transformation 

(one per star tracker) 

d. ( TEl) = ecliptic to mean of 1950 transformation 

e. Q = navigation base to star 'tracker misalinement quaternion 

(one per star tracker) 

th 

f. STAR ^ = unit line of sight vector of the N*^ Shuttle navigation star in 

mean-of-1950 coordinates 

g. T ARGET = relative rendezvous target line of sight vector in mean-of-1950 

coordinates 

h. GAMMA = mean longitude of the Sun relative to the mean equinox of date. 

If simulated ■ sensor outputs for a star are desired, the first step is to adjust 
the line of sight of the star to simulate the effect of stellar aberration (fig. 
5-3, block 1). These phenomena result in a shift in the observed line of sight 
of a star due to the relative difference between the velocity of the observer 
and the velocity of the light from the star. The apparent shift is in the direc- 
tion of the velocity of the observer. If the rendezvous vehicle is being 
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tracked this step need not be performed. The following computations are used to 
simulate the effect of stellar aberration. 
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Figure 5-3.- Star tracker simulation functional block diagram. 
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a. Compute a unit vector parallel to the Earth's velocity vector in ecliptic 

coordinates. 

/-sin (GAMMA) \ 
u =( cos (GAMMA) I 

GAMMA = mean longitude of Sun at current environment time 

b. Compute the aberration correction vector 

cor = { V + s(tei] u } /c 

(tei) = mean of 1950 to ecliptic transformation 
V = vehicle velocity vector in mean of 1950 coordinates 
S = Earth's orbital speed 
G = speed of light 

0 . Add the correction vector to the star line of sight 
Sg = UNIT ( STAR j,j + cor ) 

Before adding error to the adjusted star line of sight or the target line of 
sight it must be transformed into the appropriate star tracker coordinate sys- 
tem. The onboard navigation software uses an ideal navigation base to star 
tracker transformation; however, the actual, or simulated, transformation will 
be in error by an amount defined by the input navigation base to star tracker 
misalinement quaternion, Q. Q is defined as 
Q = (Qq, Q) 

Qq = cos co/2 
Isl = sin ti)/2 

unit (Q) = direction of misalinement error 
w = magnitude of misalinement error 
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The two sources of misalinement error, w, for the star tracker are the naviga- 
tion base to star tracker misalinement and internal star, tracker bias . , The 
value to be used for ti), therefore, is the RSS of the magnitudes of these two 
errors. The actual navigation base to star tracker transformation is computed 
in block 2 of figure 5.3 by the following equation. 

1-2(02^ + Q32) 2(QiQ2 - Q0Q3) 2{QiQ3 + Q0Q2) 

I^TNsJactual = 2^^201 + Q0Q3) 1-2(Q3^ +. Ql^) 2(0303 - QoQi) ^TNsJi^eal 

2(030]^ - Q0O2) 2(0303 + 0 qQ]_) 1 - 2 (Q]_^ + 03^) 

(TNS]^^g^-j^ = ideal navigation base to star tracker transformation 

(one per star tracker) 

[TNS)gg.(.yg 2 = actual simulated navigation base to star tracker transformation 
(one per star tracker) 

For a single mission simulation the navigation base to star tracker misaline- 
ment, 0, will be constant; therefore, the transformation t^^^^actual 
only be computed once. 

The apparent star or target line of sight can now be rotated into the appropri- 
ate star tracker coordinated system (block 3). 

(rBB](TIB]s^ 

(tbn) = vehicle body to navigation base transformation 
[tIb) = mean of 1950 to vehicle body transformation. 

The error model (block 4 and block 5) computes and adds errors to the star line 
of sight yielding a measured star line of sight in star tracker coordinates, 
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The following equations are then used to compute the equivalent star tracker 
output angles for the measured line of sight (block 6) . 

H = tan"^S^2/Sm3) 

V = tan"^(-Sin2/Sjj,3) 

These equations are derived from the relation between the star tracker measured 
output angles, H and V, and the three star tracker coordinate system axes, h, 


5-16 



ORIGINAL PAQE IS 
OF POOR quality 


BORES IGHT AXIS 

0 




DEFINED AS +h, +v, BORESIGHT 


h, V STAR TRACKER AXIS 

H, V OUTPUT DEFLECTION ANGLES 




Figure 5-4.- Star tracker coordinate system. 
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5 . 3 ERROR MODEL 

The star tracker error model computes simulated measured star tracker output 
angles by adding errors to the line of sight of the object being tracked. This 
error model is also designed to be used in conjunction with the star tracker 
data accumulation software (sec. 5.1.2). 

The sources of star tracker error can be divided into. two groups. The first 
group contains errors that will be referred to as the random errors. These 
errors appear randomly distributed when sampled at a high frequency or at the 
star tracker software data sampling rate of 6.25 Hz. The sources of these 
errors include: 

a. Noise 

b. Digital resolution- 

c. Variations in input voltage 

The remaining star tracker errors remain fairly constant over the period of time 
required for the star tracker software to accumulate 21 sightings on 1 particu- 
lar star. These errors, although random, change noticeably only over long pe- 
riods of time and for different stars. For a set of 21 sightings on 1 star for 
a period of 3*36 seconds these errors would not change significantly between 
each individual sighting. The sources of these bias errors include; 

a. Variations in star intensity 

b. Variations in temperature 

c. Earth's magnetic field 

d . Time lag 

e. Mechanical stability 
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f. Lens stability 

g. Calibration 

h. Long-term drift of electronic components 

This model adds just two errors to each star sighting: one being the ESS of all 
the random errors and the other the RSS of all the bias errors. The random 
error is recomputed and changed for each individual star sighting*, while the 
bias is computed once for each set of 21 star sightings and added to each indi- 
vidual sighting of that set. 

This model adds the errors to a line-of-sight vector by rotating it in a paral- 
lel plane by an amount equal to the error. The amount of rotation is computed 
from a normally distributed random number, the direction of rotation is deter- 
mined by a uniformly distributed random number, and the actual rotation is 
performed by vector quaternion multiplication. This method of adding on the 
error assumes that the random error is both Gaussian in distribution and 
Isotropic. 

The following is a detailed description of the computations used in blocks 4 and 
5 of figure 5-3. 

a. If this is the first sighting of the 21 sightings, compute the bias error 
and construct the bias error quaternion for this set of sightings; 
otherwise, proceed to step 2. 

(1) Construct the bias quaternion so that it is perpendicular to the ac- 
tual line of sight and parallel to the Y-Z plane. 

V = UNIT X i) 


5-19 



79FM5 


(2) Compute the direction of the bias quaternion and rotate it about the 
actual line of sight by an amount equal to the computed direction. 

R = (Rq, R) 

Rq = cos Cir/2 
R = (sin CTT/2) 

C = random uniformly distributed number 

A value for C is supplied by a uniform distribution random number 
generator such that 0 ^ C < 1; therefore, 0 < CiT < ir. 

V = R* V R 

(3) Compute the magnitude of the bias quaternion. 

QB = (QBq, 

QBjj = cos (E * BIAS/2) 

QB = {sin (E ♦ BIAS/2) } V 

BIAS = RSS of the star tracker la bias errors 
E = random normally distributed number 

The value of E is supplied by a normal distribution random number 
generator such that -3 £ E < +3; therefore, 

-3 * BIAS < E * BIAS < +3 * BIAS. 

Add the bias for this set of sightings to the actual line of sight (block 
4, fig. 5-3). 

Sb = QB* S^^ QB 

Compute the random error quaternion for this individual sighting. 

(T) Construct the random error quaternion so that it is perpendicular to 
the vector and parallel to the X-Y plane. 
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V = UNIT (Sjj X i) 

(2) Compute the direction of the random error quaternion and rotate it 

about the line of sight by an ^ount equal to the computed direc- 

tion . 

R = (Rq, R) 

Rq = cos CiT/2 
R = (sin cn/2) 

C = random uniformly distributed number 
■ V = R* Si3 R 

(3) Compute the magnitude of and form the random error quaternion- 
QE = (QEq, 

QEq = cos (E * ERROR/2) 

QE = {sin (E * ERROR/2)} V 

ERROR = RSS of the la star tracker random errors 
E = random normally distributed number 

3. Add the random error for this sighting (block 5). 

Sm = QE * Sb QE 

rhe vector points in the simulated star tracker sensed direction of the 
star or target being tracked. This line-of-sight vector is converted into the 
equivalent star tracker output angles in block 6 using the appropriate equations 
Dresented in section 5-2. 
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5.4 ERROR SOURCES AND CONSTRAINT VALUES 

The star tracker bias error is listed in table 5-1. This error is modeled by 
the nonvarying star tracker to navigation base misalineraent quaternion. The 
star tracker measurement random error sources and the measurement bias sources 
are also listed in table 5-1. Table 5-II gives representative values for 
I^TNsjideal’ which are the I_LOAD navigation base to star tracker transformations. 
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TABLE 5-1.- STAR TRACKER MISALINEMENT AND STAR SIGHTING ERRORS 


Star tracker alineraent uncertainty (Icr) 

Star tracker to navigation base misali-nement 

“y 

"z 

star sighting errors (2-axis, ia) 

Random star position errors (ERROR) 

Star sighting biases (BIAS) 


22.8 sec 
22.8 s^ 
20.6 sec • 


15.0 sec 

60.0 sec 
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TABLE 5-II.- ‘IDEAL. NAVIGATION BASE TG STAR TRACKER TRANSFORMATIONS' 


base to -Z star 

tracker transfo”'"='^^*'^’‘' 

Component 

Values 

1.1 

-0,0056491 

1.2 

.9994101 

1,3 

-.0338744 

2,1 

.9894338 

2,2 

.0006786 

2,3 

-.1449833 

3,1 

-.1448747 

3,2 

-.0343355 

3,3 

-.9888540 

Ideal navigation base to -Y star 

tracker transformatic 

Component 

Values 

1,1 

-0.9662658 

1,2 

-.1833851 

1,3 

.1808317 

2,1 

-.1839513 

2,2 

.0000000 

2,3 

-.9829353 

3,1 

.1802558 

3,2 

-.9830411 

3,3 

-.0337339 
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6.0 CREW OPTICAL ALINEMENT SIGHT 

6.1 GENERAL DESCRIPTION 

The Shuttle crew optical alineraent sight (COAS) is an optical sighting instru- 
ment that will be used to provide both backup. IMU platform alinement and 
rendezvous vehicle tracking capability. Since the operation of the Shuttle star 
trackers requires a reasonably well known IMU reference, the COAS will be used 
when the orientation of all three IMU platforms is. unknown. The COAS must also 
be used in the event of double star tracker failure. The COAS instrument can be 
mounted in two different locations. When mounted at the forward windows its 
line of sight is in the +X body axis direction, and when installed at the 
rendezvous window it points in the direction of the Shuttle -Z body axis. 

To use the COAS for the purpose of IMU alinement, the crewman must first exeute 
the COAS sighting software via a keyboard entry item. The crewman also inputs 
the identification number of the star to be sighted and the identity of the COAS 
instrument being used. He then maneuvers- the vehicle so that the selected star 
is centered in the reticle of the appropriate COAS instrument. When the crewman 
is satisfied that the star is sufficiently centered he takes a mark by 
depressing a hardware mark button (ADI attitude reference button). When 
the mark is taken, the onboard Shuttle software computes a line-of-sight 
vector to the star in IMU platform coordinates based on the calibrated COAS 
pointing direction and IMU attitude information. After sighting data have 
been taken for two stars the orientation of. the IMU platform with respect 
to the M50 coordinate system can be determined. 

Rendezvous vehicle tracking using the COAS is performed in a similar manner. In 
this case , however , a pair of COAS calibration angles and vehicle body attitude 
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information at mark time are provided for incorporation into the onboard 
rendezvous navigation filter for use in relative state estimation. 

Before the COAS instrument can be used for star or target sighting, its pointing 
direction in both navigation base and body coordinates must be accurately deter- 
mined or calibrated. This must be done while the IMU is accurately alined. The 
crewman executes the COAS calibration software and inputs the COAS instrument 
identity and the identity of the star to be sighted. When the mark is taken on 
the selected star the onboard Shuttle software computes a line of sight for the 
COAS in both navigation base and body coordinates, based on the known IMU 
•reference and the known position of the star. The COAS calibration will 
nominally be performed for each COAS instrument as soon as is practical after 
orbital insertion. This ensures that the COAS is calibrated before any possible 
IMU or star tracker failures. Furthermore, the quality of the IMU alinement is 
maximum at this time. 
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6.2 SIMULATION APPROACH 

The functional block diagram in figure 6~1 illustrates the method used in the 
simulation of the COAS instrument . The diagram describes not only the functions 
performed by the COAS error model but also those functions performed by an IMU 
model and onboard flight software. (Note that the COAS model alone- is not suffi- 
cient for generating simulated COAS measurements since the model merely defines 
the body attitude at mark time.) In the following subsections the COAS model 
and the onboard software will be discussed separately. The IMU model is 
described in section 2.0. 

6.2.1 COAS Model 

a. S, the unit line-of-sight vector of the star to be sighted in M50 
coordinates 

b. T, the unit line-of-sight vector of the rendezvous target in M50 
coordinates 

c. C, the actual COAS unit line-of-sight vector in the vehicle body coordi- 
nates (one per COAS mounting location) 

d. (tei], the eoliptic-to-M50 transformation matrix at the time of the 
sighting mark 

e. GAMMA, the longitude of the Sun at the sighting mark time 

If simulated measurements for a star are desired, the first step is to adjust 
the line of sight of the star to simulate the effect of stellar aberration (fig. 
6-1, block 1). These phenomena result in a shift in the observed line of sight 
of a star due to the relative differences between the velocity of the observer 




Figure 6-1.“ COAS sighting functional block, diagram. 
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and the velocity of the light from the star. The apparent shift is in the direc- 
tion of the velocity of the observer. If the rendezvous vehicle is being 
sighted this step need not be performed. 

The equations necessary for computing and adding the aberration correction are 
presented in section 5.2. The error model (block 2) then computes and adds 
error to the star or target line of sight yielding a measured star line of 
sight in M50 coordinates, 

In order to generate IMU gimbal angles at the mark time the vehicle body orienta- 
tion must be determined at that time (block 3)* This involves a vehicle maneu- 
ver to aline the COAS instrument with the measured star line of sight. The ma- 
neuver can be simulated as a body rotation about a single axis until the 
measured star line of sight is colinear with the actual COAS pointing direction 
vector. The two line-of-sight vectors are first rotated into a common coordi- 
nate frame . 

Sb = [tie) ^ 

Sb = measured star line of sight in body coordinates 
(tie) = M50 to vehicle body transformation matrix at mark time 
^ = [TBN)^Ca 

Cb = actual COAS pointing direction vector in body coordinates 
(ten] = actual vehicle body to navigation base transformation 
C^ = actual COAS pointing direction vector in navigation base coordi- 
nates computed from the ideal COAS line of sight by the error 
model (block 6) 
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The axis of rotation is perpendicular to both Sg and Cg. 

E = UNIT (Cg X Sg) 

R = unit vector parallel to the axis of rotation 

The M50 to vehicle body transformation at the sighting mark time is, therefore, 

(tie) = (Cr • R ’ (C„ X R))(S„ -.R • (S„ X R))'^(tIb} 

With the M50 to vehicle body transformation matrix at mark time, IMU gimbal 
angles can be generated by using the IMU model (sec. 2.0, block 4). 

The calibrated COAS line of sight and the COAS calibration angles must also be 
generated for Shuttle software use (block 7, fig* 6-1). The computations 
performed in block 7 need only be executed once per single mission simulation to 
model the COAS calibration error. The simplest and most straightforward method 
for modeling the calibration error is to actually simulate the COAS calibration 
process . This is accomplished by using the COAS model and the onboard calibra- 
tion software (fig. 6-2). 

6.2.2 Onboard COAS Software 

The previous section described how the COAS model produces a definition of body 
attitude at mark time and the simulated COAS calibrated line of sight. The body 
attitude information, in conjunction with the IMU model (sec. 2.0), can be used 
to generate simulated measured IMU gimbal angles. This section will describe 
how the IMU and calibration data are used by the onboard software to produce 
measurement parameters that can be used by either IMU alinement or rendezvous 
navigation simulations. 
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Simulated IMU gimbal angles are used to compute the navigation base to IMU plat- 
fonn transformation at the mark time. The direction of the COAS line of sight 
in navigation base coordinates will have been previously determined for each 
COAS location; therefore, a line-of-sight vector to the star in IMU coordinates 
can be computed (block 8, fig. 6-1). 
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Sc = [TNBPcJCjjj 

Sc = measured star line-of-sight vector in IMU platform coordinates 
[tNBPC) = navigation base to IMU cluster transformation matrix at mark 
time 

Cjjj = calibrated COAS line-of-sight vector in navigation base coordinates 

The navigation base to IMU cluster transformation matrix at mark time is 
computed by the equation for [tnBPc] (sec. 5.1.3-a, block 5, fig. 6-1); however, 
in this case the IMU gimbal angles at mark time are used rather than averaged 
girabal angles . 

In order to simulate rendezvous sighting data the gimbal angle information is 
used to compute an equivalent mean of 1950 to vehicle body transformation (block 
9, fig. 6-1). 

(tMSOBd] = (tNBBD){(tCM50] [ TNBPC 

[ TNBPC ) = Navigation base to selected IMU cluster transformation matrix 

at mark time 

(tCM 50] = IMU cluster to mean of 1950 reference transformation (.constant) 

[tnBBDJ = Navigation base to vehicle body coordinate transformation 
(constant) 

(tmSOBdJ = Mean of 1950 to vehicle body transformation. 

The navigation base to selected IMU cluster transformation matrix at mark time 
is computed by the equation for (tNBPc) (sec. 5.1.3-a) using the IMU gimbal 
angles at mark time instead of averaged gimbal angles. The mean of 1950 to vehi 
cle body transformation matrix, a pair of COAS calibration angles in body 
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coordinates, and the COAS identity can then be provided for incorporation into 
the rendezvous navigation filter for use in relative state estimation. 

I 

One final software task is the computation of the aforementioned COAS calibra- 
tion angles and the calibrated COAS unit line-of-sight vector. The calibrated 
COAS unit line-of-sight vector in navigation base coordinates is computed in 
block 10 of figure 6-2. 

= (tnbpc)^ (tCM50)^ Sq ■ 

Sq = selected star unit line-of-sight vector in M50 coordinates 
(tCM5o) = IMU' cluster to M50 reference transformation matrix 
(tNBPC) = navigation base to IMU transformation matrix 
Cjjj = calibrated COAS line of sight In navigation base coordinates 

The navigation base to selected IMU cluster transformation matrix at mark time 
is computed by the equation for (tNBPC] (sec. 5.1.3-a) using the IMU gimbal 
angles at mark time in place of averaged gimbal angles. COAS calibration angles 
in vehicle body coordinates are also computed for the rendezvous navigation fil- 
ter (block 10. fig. 6 - 2 ). 

Ce = (tNb] Cjj, 

(tNb) = navigation base to vehicle body transformation 
a_z = -tan“^ 

= tan“^ (Cei/Ce3) 

= tan'^ (Cgg/Cgi) 

6+x = 
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6.3 ERROR MODEL 
6«3'1 Sighting Error 

The COAS instrument sighting error model generates a, simulated measured line- 
of~ sight vector by adding error to the apparent star or target line of sight.. 
This model adds error to a line of sight by rotating it in a parallel plane by 
an amount equal to the error. The amount of rotation is computed from a nor- 
mally distributed random number, the direction of rotation is determined by a 
uniformly distributed random number, and the actual rotation is performed by 
vector quaternion multiplication. This method of adding on the error assumes 
that the random error is both Gaussian in distribution and isotropic. 

The following is a detailed description of the computations used in block 2 of 
figure 6-1. 

a. Compute the random error quaternion for this sighting. 

( 1 ) Construct the error quaternion so that it is perpendicular to the 
vector ^ and parallel to the X - Y plane. 

V = UNIT (S^ X i) 

(2) Compute the direction of the error quaternion and rotate it about 

the line of sight by an amount equal to the computed direction. 

R = (So, R) 

Rq = cos dTT/2 
R = (sin dTT/2) 

d = random uniformly distributed number 
A value for d is supplied by a uniform distribution random number 
generator such that 0<d<l; therefore, 0<dTr<7f. 
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V = R* V R 

(3) Compute the magnitude of and form the error quaternion. 

Q = (Qq, Q) 

Qq = cos (e * ERROR 1^/2) 

Q = {sin(e*ERR0Ro/2)} V 

ERRORq = ICT COAS instrument random error 

e = random normally distributed number 

The number, e, .is supplied by a normal distribution random number 
generator such that -3 < e < +3; therefore, 

-3*ERRORq < e » ERROR Q < +3*ERR0 Rq. 

b. Add the random error for this sighting 

Sm = 5* « 

The vector Sjjj points in the simulated measured direction of the star or target 
being sighted. 

6.3*2 Structural Instability Error 

After the COAS has been calibrated, its line of sight in body coordinates will 
change due to random structural instabilities. This COAS mounting uncertainty 
is modeled using the same method that is used for adding the sighting error to 
a star line of sight. The following is a brief description of the computations 
performed in block 6 of figure 6-1. 
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a. Compute the COAS mounting uncertainty error quaternion. 

(1) Construct the error quaternion so that it is perpendicular to the 
vector C and parallel to the body X-Y plane. 

V = UNIT (C X i) 

(2) Compute the direction of the error quaternion and rotate it about the 
ideal COAS line of sight, C, by an amount equal to the computed direc- 
tion , 

R = (Rq, R) 

Rq = cos d'lT/Z 
R = (sin dTT/2) C 

d = random uniformly distributed number 

V = R* V R 

(3) Compute the magnitude of and form the error quaternion. 

Q = = (Qo» Q) 

Qq = cos (e * ERROR^/2) 

Q = sin (e * ERROR .|/2) V 

ERROR ^ = 1 a COAS mounting uncertainty due to random structural 
instabilities 

e = random normally distributed number 

b. Add the error to the ideal COAS line of sight vector. 

^ = Q* C Q 

The vector C^ defines the actual pointing directions of the COAS instrument in 
vehicle body coordinates. 
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6.4 ERROR SOURCES AND CONSTRAINT VALUES • 

The sources of COAS , measurement error are listed below. 

a. COAS sighting error 

b. COAS calibration error 

c. COAS mounting uncertainty 

d. IMU error 

e. IMU mounting uncertainty 

The IMU errors are a subject of the IMU model description presented in section 
2.0. Of the remaining three, two of them are modeled by the algorithms 
described in section 6 . 3 . These two sources, the sighting error and the mount- 
ing uncertainty, have the following 1 sigma values (ref. 6-3). 

COAS sighting error = ERROR^ = 110 s^ ^ 

COAS mounting uncertainty = ERROR^ = 76 

As stated previously the calibration error can be modeled by actually performing 
a simulated COAS calibration (note that the calibration error is, therefore, a 
function of all of the following error sources) : 

a. (DAS sighting error 

b. Random IMU error 

c. IMU error due to drift 

d. IMU mounting misalinement error 


6-14 



79FM5 


6.5 REFERENCES 

6-1 Corson, R. W. : Shuttle Star Tracker Subsystem Operating Program Crew Opti- 

cal Alignment Sight (COAS) Simulation. JSC IN 76-FM-41, June 18, 1976. 

6-2 Star Tracker Subsystem Operating Program. Rockwell International (Downey) 
July 1, 1977. 

6-3 Holloway, T. W. : On-Orbit Flight Techniques Meeting #27 Minutes. NASA Mem- 

orandum CA5-78-48, August 31 > 1978. 


6-15 



79FM5 


7.0 RENDEZVOUS RADAR 

7.1 GENERAL DESCRIPTION 

r 

The Ku-band radar /communication system is a dual-purpose Orbiter subsystem that 
is used during orbital operations. It functions as either a wideband communica- 
tions system for data interchange with the ground or as a rendezvous radar to de- 
tect and track detached payloads or other spacecraft before or during rendezvous 
maneuvers. The functions of communications and radar are not available 
simultaneously . 

The rendezvous radar provides automatic target detection, acquisition, and track- 
ing. It is a target location and parameter estimation device that assists the 
navigation filter in reducing" target location uncertainties during rendezvous 
maneuvers. The targets are Earth- orbiting satellites that are stabilized in 
three axes. The radar tracks a target automatically and supplies range, range 
rate, radar shaft (roll) and trunnion (pitch) angles, and Inertial angle rates 
to both the navigation filter and the crew displays. All output parameters are 
provided simultaneously during the tracking phase. External assistance and 
operator intervention are normally required only for functions that initiate the 
radar search operation. 

Upon the initiation of search, the radar performs automatic search about a direc- 
tion (vector) furnished to the radar by the GPC and indicates that search is in 
progress for a maximum of 1 minute. If target detection occurs, the search 
program stops, the search Indication is removed, and the radar indicates that de- 
tection has occurred. 
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The radar operates in two main modes: passive (no transponder on target vehicle) 
and cooperative (transponder operating). The range of the unit is enhanced by 

t 

a factor of 25 when in the cooperative mode, as well as increased measurement 
limits, while the error characteristics remain essentially the same. 
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7.2 SIMULATION APPROACH 

The functional block diagram in figure 7-1 illustrates the method used in the 
simulation of the rendezvous radar hardware. It is assumed that the following 
five inputs are available from an environment simulation. 

a. ^3 = Shuttle Orbiter position vector in M50 coordinates 

b. V = Shuttle Orbiter velocity vector in M50 coordinates 

o 

c. R„ = rendezvous vehicle position vector in M50 coordinates 

R 

d. V = rendezvous vehicle velocity vector in M50 coordinates 

R 

e. [tIb] = M50 to Shuttle vehicle body transformation matrix 

The following discussion summarizes the computations performed in each of the 
blocks in the radar model functional block diagram (fig. 7-1). 

a. Block 1 - compute the Shuttle relative position and velocity vectors of the 
rendezvous vehicle in the mean of 1950 coordinate system. 


R 





^ = rendezvous vehicle position vector relative to the Shuttle in M50 coordi- 
nates 

V = Shuttle-relative rendezvous vehicle velocity vector in M50 coordinates 


b. Block 2 - compute the mean of 1950 to radar sensor coordinate transforma- 
tion, (tisJ. 

(tis) = (tbs) (tib) 

(tbs) = Shuttle body to radar sensor coordinate transformation 
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Figure 7-1.- Rendezvous radar model functional block diagram. 
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The rendezvous radar coordinate system is^ oriented with respect to the Shut 
tie body as illustrated in figure 7-2, and the body to radar sensor transfo 
mation matrix is defined as 

0 0-1 

0 10 

10 0 

c. Block 3 - compute the target range, r. - 

r = |R| 

d. Block 4 - compute the target range rate, r. The range rate is the compo- 
nent of V that is parallel to R; therefore, 

r = (R * V)/r 

Extract the. directions of the radar sensor system coordinate axis, X, Y, 
and Z, from the M50 to sensor transformation. 

[ X . , Y . Z } = (tIS]'^ 

• ^ • 

Compute the component of the range vector, R, that is perpendicular 
to the radar shaft axis (fig. 7*3). 

P = R - Z (Z • R) 


[tbs) = 
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e. Block 5 - compute the shaft angle, 'F, and the trunnion angle, 0. The tan- 
gent of the shaft angle is simply the component of the range vector parallel 
to the sensor Y-axis divided by the component of the range vector parallel 
to the sensor X-axis; therefore, 

= tan"’’ { (R • Y)/(R ’ X)^ 

Similarly, the tangent of the trunnion angle is the component of the range 
vector parallel to the sensor Z-axis divided by the component of the range 
vector parallel to the perpendicular, P. 

Q = tan"'' { |p1 (R • Z)/(R * P) } 

» 

f. Block 6 - compute the inertial roll rate, and the inertial pitch rate, 

0. The roll rate is equivalent to the component of the relative velocity 
vector perpendicular to both R and Z divided by the length of the perpen- 
dicular P. 

^ • (Z X R ) }/ (P • P) 

Similarly, the pitch rate is equal to the component of the velocity vector 
that is perpendicular to both the range vector and the trunnion axis divided 
by the magnitucie of the range; therefore, 

0 = y • {R‘x (Z X R)>/(R * R |p|) 
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g- Block 7 - Rendezvous radar error model. Once the measurements have been 
computed they must be tested to ensure that they are realistic with respect 
to the limits of the instrument. If any measurement is bad (i.e., violates 
the constraints listed in table 7-1), then all measurements are bad and the 
output data good flag is set to low or zero. This occurs with the 
rendezvous radar since all measurements depend on having lock on the target. 
The ideal radar output parameters are then corrupted using the equations in 
section 7.3* 
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7.3 error model 

The computation of errors is' 'divided into two parts. .At the beginning of each 
phase (and only then) the random bias for each measurement type is computed. 1 
form of this computation is 

RANDOM NUMBER *' 1-SIGMA BIAS = RANDOM BIAS 

where the random number is Gaussian with zero mean and standard deviation of 1, 

The noise is computed every pass through the model, and follows the same for- 
mula: 

RANDOM NUMBER * 1-SIGMA NOISE = RANDOM NOISE 

The 1-sigma values of noise and bias are selected and/or computed according to 
the requirements shown in table 7-H. 

The data shown are the current requirements. All 1-sigma values are constant 
over the measurement interval except range noise, which is range dependent. 

Depending on whether the range is greater or less than a predetermined ’’gate*' 
value, the range noise is either a percentage of the range or a fixed value, 
vdiichever is greater. 

The noise and bias are added to each ideal rendezvous radar output to simulate 
the corrupted measurement. 

^ i oi ^i 
dm = q + p + U 

i . 4-T- 

dm = measured radar output parameter 
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qi 2 j^th actual radar parameter 

B = bias for radar measurement parameter 

tl = computed noise on the radar measurement parameter 
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7 A ERROR SOURCES AND CONSTRAINT VALUES 

The radar parameter measurement limit requirements for the Orbiter rendezvous 

* • * 1 • ♦ ^ 

' . * 

radar are listed in table 7-1 . These limits are functions of both hardware and 
physical limitations.. The maximum ..range is a function of energy output and thus 
is constrained by the size, weight, and power requirements of the unit. The min- 
imum range is due to secondary reflections caused by echoes from the primary ve- 
hicle. The angle and angle rate constraints are due to limits on gimbal travel, 

travel rates during tracking, and time constraints placed on the acquisition of 

/ 

the target during the acquisition mode. (It should be noted that the angle rate 
constraint is tested by components, not by root sum squared, since the azimuth 
and elevation angle gimbals are independent.) Range rate is constrained by the 
pulsing frequency of the radar (rate at idiich it receives and processes 
information) . 

The radar parameter measurement error requirements are presented in table 7-H • 
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TABLE 7-1.- RADAE PARAMETER MEASUREMENT LIMITS 


Parameter 

Passive mode 

Active mode 

LOS range 

100 ft (30 m) to 
10 n. mi. (19 km) 

100 ft (30 m) to 
300 n. mi. (560 km) 

LOS range rate 

W8 ft/sec 
(45 m/sec) 
closing, 75 ft/sec 
(23 m/sec) opening 

1500 ft/sec 
(457 m/sec) 
closing, 300 ft/sec 
(91 m/sec) opening® 

LOS pitch angle 
(relative to 
Orb iter minus-Z 
axis)*^ 

Plus or minus 
30 deg 

Plus or minus 
30 deg 

LOS roll angle 
(relative to 
Orbiter minus-Z 
axis)^ 

Same as pitch 

Same as pitch 

Pitch Inertial 
angle rate 

Plus or minus 
20 mr/sec 

Plus or minus 
20 mr/sec 

Roll inertial 
angle 

Same as pitch 

Same as pitch 


^The LOS range rate limits versus range shall be as follows: : 

1500 ft/sec closing, 300 ft/sec opening at 300 n. mi. 

700 ft/sec closing, 300 ft/sec opening at 200 n. mi. 

300 ft/sec closing, 300 ft/sec opening at 100 n. mi. 

200 ft/sec closing, 300 ft/sec opening at 50 n. mi. 

148 ft/sec closing, 75 ft/sec opening at 10 n. mi. to 100 ft. 

^Angle measurements capability shall be maintained over the full coverage ca- 
pability of the antenna with accuracy degraded to levels no less than the 
equivalent of communication angle measurement accuracy. 
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TABLE 7-II.- radar PARAMETER , ALLOWABLE MEASUREMENT ERRORS 


Param^eter 


1 


Random error 

(3a) 


Bias error 

(3a) 


LOS range 


LOS range rate 


LOS pitch angle 


1 ft/sec (0.3 m/sec) +1 ft/sec (0.3 m/sec) 
or 1 percent of range 
rate, whichever is . 
greater 

8 mrad +3' deg*^ 

(0.458 deg) 


LOS roll angle 


8 mrad 
(0.458 deg) 


+3 deg^^ 


Pitch inertial 
angle rate^ 

Roll inertial 
angle ° | 

I 


0.14 mrad/ sec 
(0.008 deg/sec) 

0.l4 mrad/sec . 
(0.008 deg/ sec) 


+0.14 mrad/sec 
(0.008 deg/sec) 

+0.14 mrad/sec 
(0.008- deg/sec) 


^See: figure 7-4. 

°Notj including target effects. 

°Seej figure 7.5. 

*^Includes 2-degree requirement value plus a structural uncertainty 
of 11 degree. 
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Figure 7-4.- Three-sigma random range error. 
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8.0 RADAR ALTIMETER 

8.1 GENERAL DESCRIPTION 

The radar altimeter (RA) is a C-band altimeter system that provides a measure- 
ment of ground-relative altitude of the antenna system of an airborne vehicle 
from 5480 feet altitude to touchdown. The RA requires no ground station infor- 
mation as is the case with the other navigation sensors. The user simply 
applies power to the RA and the system provides altitude from the antenna to the 
nearest point on the terrain. Thus, the RA reads 4+0.5 feet when the 
Orbiter's wheels are on the ground. 

To obtain altitude information the airborne transmitter illuminates the terrain 
below the aircraft with periodic 4.3-GHz pulses. The airborne receiver then mea- 
sures the time delay between transmission and the first echo received (denoted 
leading edge detection). The two-way propagation time is then multiplied by 
one-half the propagation velocity to compute the aircraft's altitude above- the 
nearest terrain point. The altitude measurement is then sent to the onboard com- 
puter at a 6.25“Hz rate for use in guidance computations and attitude vertical 
velocity indicator (AVVI) processing. A repetition rate of 8.5 kHz for the 
transmit-receiver sequence allows the RA to track altitude change rates of up to 
2000 ft/sec. 

The radar altimeter will operate accurately over any terrain type that exhibits 
a coefficient of reflection of 0.05 or greater. The RA may, therefore, be used 
over water or ice as well as land. The RA may be used in heavy rainfall of up 
to 40 mm/hour. 



8-1 



79FM5 


8.2 SIMULATION APPROACH 

The Shuttle radar altimeters lock on and yield valid altitude measurements at an 
altiude of 5^80 feet. Altimeter data are not utilized, however, until the 
sub-Orbiter point is within a specified distance from the runway threshold. 

This limit, denoted R_G0_RAMP, may be found in the navigation I-load (ref. 8-6). 
The current value for R_G0_RAMP is 4000 feet. Assuming a 1.5~degree glide- 
slope after final flare, and an aim point 1500 feet beyond the runway threshold, 
radar altimeter data will first be incorporated in the nav. state vector sent 
to guidance at a nominal altitude of 144 feet (= (1500+4000) sin 1.5 deg). 

After RA lockon, altitude measurements are simulated via the altitude-dependent 
error model (sec. 8.3). 


8-2 



79FM5 


8.3 ERROR MODEL 

Radar altimeters are, in general, subject to the following four potential error 
sources : ^ 

a. Averaging of the terrain illuminated by the EM beam 

b. Multipath effects 

e. Error in pointing of the altimeter EM beam or terrain slope 
d. Errors inherent in the altimeter measurement equipment 

The technique of leading edge detection (utilized by the Space Shuttle radar al- 
timeters) effectively eliminates the second' error source. Leading edge detec- 
tion' yields altitude above the highest or nearest terrain point. The first and 
third error souces (terrain averaging and terrain slope) are not significant 
near airstrips where the Shuttle is likely to land. In special cases when the 
Orbiter makes its final approach over uneven terrain (as is the case for certain 
approaches to Vandenberg) it is likely that the threshold point at which RA mea- 
surements are incorporated will be moved closer to the runway threshold so that 
altimeter data over rough terrain will not be used. 

Measurement equipment error constitutes the most significant error in the Shut- 
tle RA system. Limited flight-test data gathered by, the manufacturer 


^A fifth RA error source has emerged during the Approach and ‘Landing Test 
(ALT) program. ALT free-flight test data have indicated erroneous altitude 
transients (from RA #1) of 45 feet when the Orbiter is below 100 feet 
altitude on final approach. To correct this problem for the orbital 
flight tests (OFT), the RA antennas will be repositioned so that Interference 
between the two RA*s and the landing gear is minimized. The error model 
presented here assumes that this RA altitude transient problem will be 
corrected for OFT. 
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irom wo aj-timeT^ers riown in a aircrart indicate that Dias accounts for ap- 

proximately 97 percent of the total RMS system error. The residual RMS system 

2 2 1/2 

error of 25 percent ( = {100 percent - 97 percent ) ) is attributed to random 

noise. Radar- altimeter error is modeled as follows: 

h^ = h + bjj U^ + CTjj 

where h„ is the corrupted altitude measurement 

h is the true altitude 'of the radar altimeter antennas above the 
local terrain 

^h’^h bias and noise figures given in section 8.4 

u^ is a unity variance, zero-mean Gaussian random number, sampled 
once per flight for each altimeter 

Uo is the same as u,, but sampled once per measurement for each 
altimeter . 

Finally, the altitude measurement is quantized to the nearest foot via: 
h = 0.5 ft + q(h^) 

where q() is the quantization operator that truncates the argument- to the least 
significant bit for the RA altitude word. 
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8.4 ERROE SOURCE VALUES ; 

The Orbiter carries two Honeywell model HG7193C1 radar altimeters: one each i 

the 1 and 2 forward avionics bays. - The following error model statistics are 
based on the Shuttle RA procurement specification (ref. 8-3)* 


Bias (IcT): Value 

0 < h < 10 ft 1.552 ft 

10 < h < 66 ft 2.522 ft 

66 < h < 5476 ft .0291 h 

Noise (Ig): 

0 < h < 10 ft 0.389 ft 

10 < h < 66 ft .632 ft 

66 < h < 5476 ft .0075 h 


Quantization 


1 ft 
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9.0 TACM 

9.1 GENERAL DESCRIPTION 

Tacan (tactical air navigation) is a military L-band air navigation system that 
provides slant range and magnetic bearing from user aircraft to the ground bea- 
con. Tacan is an internationally accepted system with over 1300 ground stations 
worldwide, of which over 800 are within the United States. 

The user selects a desired ground beacon by tuning to one of the 126 channels in 
the frequency band from 1025 MHz to 1150 MHz (ref. 9-1). At intervals of 30 sec- 
onds each ground beacon broadcasts its unique station ID in audible Morse code 
tones to allow positive station identification by the crew of the user aircraft. 

To obtain slant range information an airborne interrogator transmits a pair of 
pulses omnidirectionally on a specific frequency appropriate to a selected 
ground station , then searches for a returning pair of pulses from the ground sta- 
tion transponder (receiver, retransmitter) on a frequency +63 MHz from the 
interrogation frequency. The airborne unit then determines slant range distance 
by measuring the two-way propagation time, subtracting 50 Dsec to account for 
ground transponder delay, and multiplying by one-half the propagation velocity. 
This one-way distance is then provided to the onboard computer in units of nauti- 
cal miles . The slant range function of the Tacan system is often called dis- 
tance measuring equipment (DME) . Each Tacan ground station can provide DME ser- 
vice to as many as 100 aircraft simultaneously. 

Bearing Information is transmitted from the ground station whether or not DME 
interrogation pulses are being received. The transmitted signal is a nine-lobe 
cardioid radiation pattern at a fixed 15 revolution/second rate. As one 
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particular lobe passes through magnetic north, a particular string of pulses 
(code) is transmitted. A different string is transmitted for each of the other 
eight lobes. Bearing is determined by the receiver by measuring the time differ' 
ence between the most recent reference pulses and a signal strength maximum on 
the rotating cardioid of known rate. The mechanization permits a 360° bearing 
sector to be divided into 9 identifiable 40° sectors, for a nine-fold increase 
in accuracy in determining bearing as compared to single lobe systems.^ 


^It has been shown during the Approach and^Landing Test program that Tacan 
bearing data are subject to transitory 40 spikes. Navigation designers 
should know that such transients may occur during OFT entry, and should 
be aware of their effects on the navigated state vector. 
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9.2 SIMULATION APPROACH 

Tacan provides a slant range measurement (p) from the Orbiter to the ground sta- 
tion, as well as a measurement of magnetic bearing (0) from the Orbiter to the 
ground station (w.r.t. magnetic north at the ground station). Signal strength 
and low elevation constraints are reflected in the range and bearing lockon 
flags (Piook» ®lock^» which are provided by the Tacan simulation. (Note: the 

user program must test for the high elevation constraint (elevation angle > 45°) . 
or the cone of confusion, in which Tacan bearing -data becomes unreliable.) This 
Tacan simulation includes the automatic antenna selection logic between upper 
and lower antennas for each navset that is found in the TACAN SOP (ref. 9-2). 

An environment simulation provides the true Orbiter position and velocity, ^Qpg 

and VQpg, as well as the Orbiter 's attitude as described by the transformation 

BODY 

from Earth-centered inertial (ECI) to body coordinates, ^ECI 

This Tacan simulation requires the following data for each ground station. 

a. Ground station position, RgND 

b. Magnetic variation at ground station, UyAR 

c. Ground station equipment type 

These data are available from the Environmental File, vdiich is prepared- and 
maintained by the Department of Defense Electromagnetic Compatibility Analysis 
Center (ECAC) in Annapolis, Maryland. Limited data concerning Tacan 's in the 
KSC and Edwards areas are available from references 9-3 and 9-4. 

It is assumed here that the Tacan tuning schedule has been chosen to avoid the 
following conditions: 
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a. Ground stations obstructed by surroundings (mountains, buildings, etc.) 

b. Cochannel interference caused by simultaneously receiving two ground sta- 
tions on the same channel (ref. 9-3, and 9-5). 

Figure 9-1 presents a block diagram for a Tacan simulation. The functions of 
each block are as follows: 

a. Block 1 - computes five parameters specifying the Orbi ter /ground station 
geometry. Equations used to derive each of the five parameters are as 
follows: 


Slant range (p) 

P = i^OEB " ^GNdI 


Line of sight vector Orb iter to ground station 


BlOS = ^IgND " ^ORB^^P 


Note: Using knowledge of the Orbiter's attitude (i.e., 

expressed in body coordinates for use in block 2. 


^CI -yLOS IS 


Orbiter elevation (EL) 

EL = 90 cos“^ (BLOS • BVjGNd) 

where Uv,GND unit vector describing geodetic vertical at the 

ground station 

Magnetic bearing (9), Orbiter to ground station 

0 = tan"K(BLOS'BE,GND)/fBLOS'BN,GND)) “ ®VAR 
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Figure 9-l*~ TACAN simulation functional block diagram* 
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where i?E,GND unit vector pointing east at the ground station 

i^N GND unit vector pointing true north at the ground station 

®VAR angle variation from true to magnetic north at the 

ground station 

Radial velocity , line-of-sight vector 

^RAD = i^ORB • -LOS I 

Tangential velocity , > ^'elative to the ground station 

^TAN = I -ORB • -TAN I 

where is the horizontal unit vector perpendicular to the line of 

sight 

b. Block 2 - computes the radio frequency (RF) link loss (Llink^ using the fol- 
lowing equations extracted from reference 9-6: 

^LINK = ^ORB + ^GND + Lq + Lp + Lp + FSPL 
« 

Each of the individual terms in this equation are computed as follows: 

.^RB* Oi'biter antenna gain 

First compute the look angles from the Orbiter antennas to the ground 
station: 

a = cos"l (i^T • HLOs) 

where Uj^^ is the antenna pointing direction given in section 9.^.1. The 
particular antenna is specified by the antenna auto-select logic (block 6) 
and the navset being used. 
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ulos J.ine of sight vector in body coordinates from block 1. 

Gqrb .then computed as: 

8 

» If a < 100°, P8(a) = S C^, an dBi 

n=0 

where the Cj^'s are given in section 9.4.2. 

* After the eight-order polynomial P8(a) has been evaluated, the onboard gain 

* is computed depending on which region the look angle lies. 

* Region 1 : (a<10°) 

* Gqrb = P8(ct) + 1 2 sin (0.4m) | 

* Region 2: (10°<a £ 40°) 

* Gqbb = P8(a) 1.5 sin (0.15m) 

* Region 3: (40®<a£l00°) 

* Gqrb = P8 (a) - A3 (a) sin (F3(a) m) , 

* where A3 (a) = 2 - 0.0333 (a-40) ■ 

* F3 (a) = 0.15 - 0,000833 (a-40) 

* Region 4: (100°<a<l80O) 

* Gqrb = G4 (a) + |a 4 (oi) sin (F4(a) m) 

* where G4 (a) = -10 -0.375 (a-100) 

* A4 (a) = 0.25 (a-100) 

* F4 (a) = 0.1 + 0.000625 (a-100) 


Change 3' 
September 1980 
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* If the user desires to simulate the effects of variations in the on-board anten- 

* na patterns, then a- random shift may be introduced in the sinusoidal scalloping 

* of the gain pattern by adding a uniformly distributed random {on [0, 2 "^)) phase 

* shift to the argument of the sine function in each region. The randcM phase 

* shift is intended to be chosen only once for each onboard antenna thereby 

* introducing a pseudorealistic pattern variation. 

GcMp: Ground station antenna gain 

GqnD = 20 log (Gg) + GtYPE 

* where Gg is the normalized antenna gain, specified by the function of 

* the elevation angle EL, as given in section 9.4.3. 

G^ype equal to 2 dB for the older ground station antenna design 
(utilized by the military) and is equal to 5 dB for the new equip- 
ment types (RTB or RTC) . Reference 9-2 gives equipment types for 
Tacan’s in the Edwards AFB area. Reference 9-3 gives similar infor- 
mation for Tacan's close to KSC. 

Lg; Cable losses 

Lg = -4.8 dB (based on the Orbiter acceptance criterion) 

Fading loss (multipath) 

If EL > 15° , Lp = 0 dB 

If EL < 15°, Lp = -(4.5 - .6 EL +.02 EL^) dB 


Change 3 
September 1 980 
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Lp : Polarization loss 

First, compute the mismatch angle (ijj) between the Orbiter and ground station 
polarization directions: 

<}> = cos-'* (Pr*ui/|PeI) 

Wll6r*6 •" ^£jOS ^ 

Hi = ^ Hlos^ ^ Hlos^I^orb ^ -losI 

Polarization loss is then computed from the mismatch angle via; 

Lp = 10 log (.5 + .49373 cos 2<}>) (dB) 

FSPL ; Free space path loss 

FSPL = 20 log (X/4 ttp) dB 

Assuming a X corresponding to f = 1.213 GHz (the shortest Tacan wavelength), 
we have: 

FSPL = -99.48 - 20 log p (dB) 

where p is in units of nautical miles 

c . Block 3 - computes received signal powers at the ground station for the 
downlink (SpQ^) and at the Orbiter for the uplink (Syp) . 

^DOWN = ^ORB + ^LINK 

where Pqrb radiated power from the Orbiter navset, assumed to be 60 

dBm (1000 watts) 


9-8 



79FM5 



9-7 


QRlGlNAt PAGE IS 

OF POOR QUALITY. 


79FM5 

Lp,: Polarization loss . 

First, compute the mismatch angle (<})) between the Oribiter and ground station 
polarization directions : 

<t> = cos“^ 

\ 

where Pp = u^rp - uj^qq cos a 

III ^^ORB * HlOS^ ^ HlOS^I^ORB ^ HlosI 

' I 

Polarization- loss is then computed from the mismatch angle via: 

? 

Lp = 10 log (.5 + .49373 cos 2'4>) (dB^ 

FSPL : Free space path loss 

FSPL = 20 log (X/4iTp) dB 

Assuming a X corresponding to f = 1.213 GHz. (the shortest Tacan wavelength), 
•we have ; 

FSPL = -99.48 - 20 log p (dB) 


where p is in units of nautical miles ‘ 
c. Block 3 - computes received signal powers at the ground station for the 

V 

downlink (Spg^) and at the Orbiter for the uplink (Sgp) . 

^DOWN = ^ORB ^linK 


where Popn is the radiated power from the Orbiter navset, assume'd to be 60 


ORB 

dBm (1000 watts) 


% 
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Llink ^ link loss computed in -block 2 

^UP = ^GND * ^LINK 

where is the power radiated by the ground station in dBm, This is 

assumed to be 6? dBm (5kW) for permanent Tacans and 58.45 dBm (700 W) mobile 
ones . 

Block 4 - computes the corrupted range measurement, p, using the range 
error model equations presented in section 9.3.1. 

Block 5 - determines range data lockon as a function of downlink received 
signal strength, Sqq^, and the elevation angle, EL, as follows: 

If EL > 10 and SpQWN > -9° <JBm, then pLocK = 

Otherwise 

plock " false 

After the range threshold criterion has been met, a period of time is 
required for the Orbiter navset to lock on to the ground station replies. 
This time increases as range rate increases, but 2 seconds may be assumed at 
range rates less than 4500 knots (ref. 9-8). 

In addition to lockon delay a memory feature exists in the hardware that 
causes the lockon indication to remain for about three-quarters of a second 
after the signal falls below threshold. 

The 1° constraint on elevation angle accounts for atmospheric effects such 
as tropospheric scatter. This is felt to be a conservative limit since the 
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hardware may actually achieve lockon below l“ elevation, but will not 
dependably do so all the time. 

The -90 dBm constraint on the received power at the ground station is due to 
the receiver threshold in the ground equipment that causes all interrogations 
with less than -90 dBm received power to be ignored. Uplink signal power 
is not a range acquisition constraint since (1) ground stations are generally 
more powerful than the Orbiter transmitters (mobile Tacan ground stations 
are less powerful, transmitting only 700 watts), and (2) airborne receivers 
are more sensitive than ground receivers, with a DME acquisition threshold 
of -94 dBm. 

f. Block 6 - selects between the two antennas (upper and lower) associated with 
each Tacan navset. 

Since range data yields more accurate position updates than bearing, it is 
desirable to switch antennas when DME is not locked-on, even though bearing 
lock mi^t be present. The Tacan Subsystem Operating Program (SOP) performs 
this antenna selection automatically, independently for each navset, in the 
following manner: 

( 1 ) The SOP will initially select the bottom antenna when a new station is 
tuned . 

(2) If, after 10 seconds (I-load), range lockon has not occurred the top 
antenna is selected. 

(3) If, after an additional 10 seconds, range lockon still has not occurred 
the bottom antenna is reselected. 
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This switching sequence continues until range lockon has occurred, at which 
time the SOP freezes to the lockon antenna. Should range look disappear for 
a 10~second period the acquisition switching sequence is resumed. 

g. Block 7 - determines bearing data lockon as a function of uplink received ' 
signal strength, Syp, and the^ elevation angle, EL, as follows: 

If EL ^ 1° and Sjjp > -90 dBm, then ®l0CK ~ TRUE 

\oCK false 

Similar to range, bearing lockon is delayed by a constant 3 seconds 
(.regardless of range rate) (ref. 9-8). Also, the data value and lockon 
indication is held for 3 seconds after loss of adequate signal. 

h. Block 8 - computes the corrupted bearing measurement, 9, using the bearing 
error model equations presented in section 9.3.2. 
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9.3 ERROR MODELS 

Tacan range and bearing errors are modeled using knowledge of the uplink and 
downlink received signal strengths, Sgp and which are computed in block 3. 

9.3.1 Range Error Model - Block 4 

Tacan range measurements are subject to biases-, bias shift, noise, and scale fac- 
tor errors as related in the following equation: 

Pq = (1 + SFE)P + bQjjj) + boRB + + ^shift + "GND * ^0 

where P^ is the corrupted range measurement 

SFE is the scale factor error 

P is the true slant range 

*^GND range bias induced by the ground station electronics 

bQRg is the range bias induced by the. Orb Iter navset electronics 

^LOC range bias induced by ground station survey errors 

^SHIFT ground station bias shift 

”GND ground station random noise 

noRB is the Orbiter navset random noise 

The range errors are modeled as independent zero-mean Gaussian random variates 
with the"^!^ values given in section 9.4.4, 

* Finally, the range is,^uantized to the nearest 0.05 nautical mile via: 

* ^*^QUANT “ I P + 0.025 0.05 (nm) 

* ■ L 0-05 

^ where INT = the greatest integer function and P = raw RNG data. 


Change, 2 
June 1980 
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9-3.2 Bearing Error Model - Block 8 

Taoan bearing measurements are subject to bias, path errors, and noise as 
described in "th 

e = e + bQfjD + bojjB + + noEB + npQUGH * ^BS 

/here 0 is the corrupted magnetic bearing measurement 

6 is the true magnetic bearing ' 

^ ^ • • 
bgjjD is the ground station bearing bias 

bQjjB is the Orbiter navset bearing bias 

”GND random bearing error induced by the ground station 

nQRB - i,® signal strength dependent random error induced, by the 
Orbiter receiver 

Trough random error due to roughness of the bearing beam 

bgg is the exponentially correlated bias due to bending and 

scalloping of the b.earing beam. The time constant (T^g) . 
of this ECRV is computed via: 

■^Bs = ^^Vrad/i^ead) + C''^tan''^tan))~^ 

?hep§ radial correlation distance 

V|j^g ±g Ih.e the line-of-sight 

Djan ia thi 

v^AN . i§ the fpethM t<b 

thg ling ©f §igh% 

igs is computed sequeafeiaiiy ia the 
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Initialization; bgs (t=0) = Ogg u 

Iteration; b^s (fc) = a'bss (t - ^t) + <^33 / 1-a^ u 

-At/Tgg 

where a = e 

u is a unit variance, zero-mean Gaussian random number 

ICf values for all the above error sources are given in section 9.4.5. 
*' Finally; the bearing is quantized to the nearest 0.1759 degree 'via; 

* BRGQojajT = INT [Q + O.O879T 0.1759 (deg) 

L 0.1759 J 

* where INT = the -greatest integer function and © = raw BRG data. 
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9.4 ERROR SOURCE VALUES 

The Orbiter carries three Hoffman ARN-84 Tacan navsets, one each in the 1, 2, 
and 3A forward avionics bays. Each receiver is fed by two L-band antennas, one 
mounted on the upper and one on the lower portion of the Orbiter 's outer struc- 
ture. Error modeling data specific to the Shuttle Tacan system follow. 

9.4.1 Orbiter Antenna Pointing Directions 


Tacan antenna 

pointing vectors 

are given 

below in body 

coordinates . 

Antenna 

Antenna 

Tacan 




number 

location 

number 

X 

I 

Z 

1 

Upper left 

1 

+0.41742 

-0.15643 

-0.89515 

2 

Lower left 

1 

+.04532 

-.49931 

+.86484 

3 

Upper center 

2 

+.42262 

0 

-.90631 

4 

Lower right 

2 

+.04532 

+.49931 

+.86484 

5 

Upper right' 

3 

+.41742 

+.15643 

-.89515 

6 

Lower center 

3 

+.03490 

0 

+.99939 


9.4.2 Orbiter Antenna Gain Pattern 

The Tacan antenna gain patterns are approximated by an eighth-order curve fit to 
the angle of offset between the ground beacon and the Orbiter antenna boresight. 
The nine coefficients specifying the curve fit are: 

Cq = -2.442764458 El 
Cl = 4.219265311 EO 
C2 = -4,327836256 E-1 
C3 = 2.344000307 E-2 
Ci| = -6^865056346 E-4 
C5 = 1.138519773 E-5 
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Cy = 5.303361563 E-10 
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9.4.3 Ground Station Antenna Gain 

Ground beacon normalized antenna gain is specified by the following function of 
elevation angle (EL) . 


» a. If EL <29 

f Gg = 0.6725 + 0.2975 sin (6) 
where 6 = 0.05265 * tt * (EL) 

« b. If 29° < EL < 90° 

* Gg = 0.212 + A2 sin (6) 

^ where A2 = 0.22 - 0.0022 (EL-22°) 

f F2 = 0.03125 + 0.0002367 (EL-290) 

^ and 6 = F2 * TV * (EL-22°) 

* If the user desires to introduce a random shift and a random amplitude variation 
^ the expression 

f V(EL) = 0.05 sin (0.1 * tt * EL + (j)) 

^ should be added to the normalized gain Gg computed above. The phase angle $ 

^ is a uniformly distributed random variable on the interval [0, 2ir) and is in- 
^ tended to be chosen only once for each ground sta" 


Change 3 
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300 ft 
100 ft 


Ground station bias shift ; (SDOWN +60) 

Noise (Id): 


Change 3 
September 1980 


9-l6a 



79FM5 


This page intentionally left blank 


Change 3 
September 1980 


9 - 1 6 b 



79FM5 


Ground station: 

®DOWN ^ “55 dBm; 20 ft 

^DOWN ^ ”55 dBms 20 — 2.33 (55 + ft 

Orbiter transceiver; 

If Sup > .-80 dBm; 42 ft 

If Syp < -80 dBm: 42-4.2 (80 + Syp) ft 

Station location uncertainty ; 33 ft (icr, per axis) 

Quantization error ; 304 ft (0.05 n. mi.) 


Scale factor error ; 100 ppm (la) 

9*4.5 Bearing Error Values 
Bias (la) 

Ground station: 0.33° 

Onboard receiver: 0.353° 


Noise (la) 

Ground station: 0.1° 

Orbiter receiver: 0.353° (Syp > -82 dBm) 

1.458° (Syp < -82 dBm) 

Quantization error ; 0.1758° (ref. 9-11) 


Path errors 

Roughness 0.3° (la) 

Bending and scalloping 0.4° (la) 

Radial correlation 10 000 ft 
distance 
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Tangential 5 000 ft 

Correlation distance 
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10.0 MSBLS 

10.1 GENERAL DESCRIPTION 

The microwave scanning beam landing system (MSBLS) provides information 
representing azimuth angle, elevation angle, and distance of user aircraft with 
respect to a ground installation located beside the runway. The ground installa- 
tion consists of two main elements; (1) an elevation station located just be- 
yond the nominal touchdown point providing good geometry for altitude measure- 
ment prior to touchdown, and (2) an azimuth/DME station located beyond the far 
end of the runway to provide data through rollout. Using knowledge of the 
ground station locations, these data are transformed into the runway coordinate 
system yielding positional information for CRT and dedicated displays, and for in- 
corporation into the Orbiter navigated state estimate. The MSBLS is a precision 
navigation sensor used- during final approach through rollout at landing sites' so 
equipped (ref. 10-1). 

The crew selects a ground installation corresponding to the desired runway by 
tuning to 1 of 10 channels via the cockpit control head. Upon reception of 
ground station solicit pulses, the onboard navsets transmit DME interrogations 
on 15.^7 GHz with a pulse pair spacing based on the channel number. Solicit 
pulses and all three information channels (Az, El, and DME) are transmitted from 
the ground stations on 1 of 10 frequency channels between 15.4 and 15.7 GHz. 

The pulse pair spacing is different for the solicits, DME replies, elevation, az- 
imuth-fly-left, and azimuth- fly-right signals facilitating navset identification 
of the data type being received. 
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DME data are decoded by -the navset using conventional transponding techniques. 

The round trip time is measured, ground station delay is subtracted, and then it 
is scaled by half the speed of light to be output in units of nautical miles. 

The ground DME antenna radiates into the entire coverage volume without scanning 
(ref. 10-1). 

The azimuth and elevation antennas are designed such that the power radiated is 
focused nearly into a plane. The azimuth antenna is oriented with this plane 
vertical and is mechanically scanned horizontally. The opposite situation 
exists for the elevation antenna. The scans result from 2.5 Hz sinusoidal oscil- 
lations, whose angular extremes are in excess of the coverage volume, but RF 
transmissions occur only during the central portion of the scan where the angu- 
lar rate is approximately constant. The sector of transmissions is plus or 
minus 15 degrees (from parallel to runway centerline) for azimuth, and zero to 
30 degrees (from the horizontal) for elevation. The solicits, DME operations, 
and angle scans are electronically s 3 mchronized for time-shared transmissions at 
200 milliseconds per complete cycle. The spacing between pulse pairs fed to the 
angle antennas is proportional to the pointing angle of the antenna at that in- 
stant. The navset tracks the amplitude modulation caused by antenna scanning 
and averages the spacing of pulse pairs received between two points of equal 
power on either side of the symmetrical antenna beam. The average spacing is 
scaled to degrees for output. Although the antenna beamwidth may be 1 degree 
(Az) and the spacing of pulse pairs is incremented only every one-eighth degree, 
angular accuracies in excess of 0.05 degree (system spec) are obtainable with 

i 

this approach. This represents less than 1 foot error per thousand feet of dis- 
tance to the ground station. 
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10.2 SIMULATION APPROACH 

Each of the triple redundant navsets aboard the Orb iter make three data words 
available to the data processing system at a 5-Hz rate. This subsystem simula- 
tion defines the algorithms necessary to model Az, El, and DME measurements 

I 

histories with respect to various error mechanisms and the behavior of the data 
validity flags. The lockon and error models utilize parameters normally avail- 
able in an environmental simulation in order to produce measurement and lockon 
histories to a relatively high degree of fidelity (subject to the accuracy to 
which these dependences are known and can be simply modeled). 

10.2.1 Block Diagram Description 

Figure 10-1 is a functional block diagram of the MSBLS subsystem simulation. 

a. Block 1 - the first task is to take the vehicle position and attitude and, 
utilizing knowledge of the Orbiter antenna and ground station locations, com- 
pute the desired line-of-sight vectors. This is done both from the ground 
stations to the Orbiter in runway coordinates and from the Orbiter to the 
ground stations in body coordinates. Slant ranges and look angles are then 
determined (sec. 10.2.2). 

b. Block 2 - the slant range from the DME site to the Orbiter and the ground 
station look angles are then checked against geometrical limits to ascertain 
whether the ground station is radiating valid data in this region (sec. 
10.3.1). This alone does not determine validity because adequate signal 
levels may not be present. 

c . Block 3 - the slant ranges and look angles are then used to compute expected 
signal levels via antenna pattern curve fits and other formulations to be 
compared to an assumed threshold for lockon determination (sec. 10.3.2). 
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Figure 10-1.- MSBLS subsystem model functional block diagram. 
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d. Block 4 - the presence of adequate RF signal and ground station coverage is 
required for each information channel validity flag. All three onboard 
navsets are assumed to behave in a parallel fashion with ^respect to data va- 
lidity (sec, 10 . 3 . 3 ). 

e. Block 5 - the look angles and signal levels are then applied to the error 
model that contains known error dependences for the ground stations, path, 
and each onboard navset (sec. 10.3.4). The three corrupted measurements for 
each data type are assumed to be representative of the performance of a non- 
failed system in the configuration experienced during the approach and land- 
ing tests (with the exceptions listed in sec. 10.3.5). 

10.2.2 Geometric Observables Computation . 

RW 

The Orbiter c.g. position' (R qq) is assumed in imits of feet in runway coordi- 
nates, and an Orbiter attitude transformation matrix (T) is assumed to convert 

a vector in runway coordinates to body coordinates by multiplication. Values to 

BODY 

be used as the Orbiter antenna location (Rqrbant^ ground station locations 
RW RW 

(^2/DME Rp;T,) are given in section 10.4,1. Lines-.-of-sight (U) from the 

ground stations to the Orbiter antenna are: 

RW r BODY , RW RW 

Uaz/dme = T" Rorbant ^ + Rcg - Raz/dme 

RW r BODY 1 RW RW 

^ Rorbant ^ + Rcg - Iel 


Lines-of-sight from the Orbiter antenna to the ground stations are then: 

BODY RW RW ^ BODY 

y. AZ/DME = ^R' AZ/DME ~ S CG ~ R ORBANT 

BODY , RW RW 1 BODY 

II. EL = ''R EL “ R CG ■ R ORBANT 
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SR AZ/DME = 

1 AZ/DME ' 

(ft) 


RW , 


SR gL = 

1 W EL 1 

(ft) 


Elevation ( 6 ) and azimuth (Y) look angles from the ground antennas to the 
Orbiter may be obtained as follows: 


Assuming: U = x + U 2 y + z 

RW RW RW 

AZ/DME = ^“'^AZ/DME 2 '^'^AZ/DME V 


(deg) 


RW . RW /iJ ^ m 

/4u 


® AZ/DME = tan [u^z/DME + («AZ/DME 2^^^ 


RW 


RW 


/rw 

^ r,T = arc tan (u^l + ^'^EL 3 ^^^ 


RW 

EL 


(deg) 


e 


RW RW RW 

EL “ tan (ug^ 3 /UgL -|) 


(deg) 


The Orbiter-to-ground station antennas look angles are similarly: 


BODY BODY BODY 

® AZ/DME = arc tan (-uj^z/dmE 3^’^AZ/DME 


(deg) 


BODY BODY /bODY BODY 

AZ/DME = arc tan (u^z/DME 2 /’^(*^AZ/DME + (^AZ/DME 3^^ ^ 
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e EL = arc tan (-UgL 3/^^^ 1 ) 
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(deg) 


BODY BODY / BODY BODY 

EL = arc tan (ugL 3^ (deg) 
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10.3 SUBSYSTEM MODELING FORMULATIONS 

This section identifies (and briefly explains) the algorithms necessary to charac- 
terize MSBLS subsystem behavior under the assumptions listed in section 10.3.5. 

The values to be used as constraints, or in detailed computation of a parameter, 
are tabulated in section 10.4. 

10.3.1 Coverage Check 

The noncorrupted, or ideal, azimuth, elevation, and DME measurements are the 
inputs to be checked against the geometrical constraint values listed in section 
10.4.2. 


TRUE 

= AZ/DME 

(ft) 

^ TRUE 

RW 

= AZ/DME 

(deg) 


RW 


TRUE 

11 

CD 

(deg) 


DME coverage is constrained by a maximum (19.6 n. mi.) slant-range value. 

If DME rj-flug < DME gp (ft) 

DME data may be valid. 

Azimuth and elevation coverage angles where data may be flagged valid are not 
the same as the previously mentioned angles where power is applied to the 
scanning antennas or the ground station specified coverage volume. 

The azimuth coverage is constrained by an azimuth limit a little in excess of 
+15°. The system sign convention is positive if the vehicle is right of 
centerline (fly left). 
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azimuth data may be valid. Similarly for elevation coverage 

If EL < EL (deg) 

elevation data may be valid. 

10.3-2 RF Lockon Cheek 

As discussed, completion of five RF links is required for acquisition of all 
three data types. They are the solicit (uplink), DME interrogation (downlink), 
DME reply (uplink), azimuth (uplink), and elevation (uplink) (ref. 10-2). A 
reduction in model complexity is possible since the solicit and DME reply links 
are identical. Furthermore, they differ from the DME downlink only by a con- 
stant (different transmission powers). The DME downlink is the weaker of the 
three such that it alone must be checked against a threshold to determine DME 
lockon. The DME downlink received signal level is also the only power 
represented in the DME error model. Azimuth and elevation powers will also be 
computed, checked, and passed to the error model. 

The RF link model approach will be to compute antenna gains, polarization loss, 
and path loss, then add a constant reflecting cable losses and transmission 
powers required for that link to get a received power. The Orbiter antenna 
gains for the azimuth and DME links are the same but a gain must be computed for 
the look angle to the elevation station. Curve fits of the four antenna gain 
contours yield gain as a function of elevation and azimuth look angles. Path 
losses are determined by the slant ranges. Polarization loss is computed by 
finding the mismatch angle between the transmitting and receiving antenna's 
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polarization sense. This angle is approximately the angle between the Orbiter Z 
axis’ and the runway X-Z plane. Received powers for the RF links may be computed 
via the formulas 


^ DME DOWN = ® T 


ORB AZ/DME DME GS 


AZ/DME 


+ G p + L pQL + L + K jjpjg (dBm) 


AZ GS ORB AZ/DME 




AZ/DME 

^ POL ^ PATH ^ AZ 


(dBm) 


EL GS ORB EL EL 

+ L poi^ + L p^jg + K gL 


^EL“^T- ■^^R 


(dBm) 


where 

S down ~ i’eceived power input to the ground receiver from the Orbiter 

(dBm) 


S AZ “ received power input to the Orbiter receiver from the azimuth 

station (dBm) 


S £L “ received power input to the Orbiter receiver from the elevation 
station (dBm) 


ORB AZ/DME 
G ^ 


Orbiter antenna gain computed as a function of the Orbiter to 
Az/DME station look angles (dBi) 


DME GS 

G p - DME ground receiving antenna gain determined by the Az/DME to 

Orbiter look angles 
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L pQj^ - polarization loss is computed by first finding the polarization 
mismatch angle (<f>pOL^* loss is then 

h POL ~ 20*L0 G.jq (cos {<j)poL^^ (dB) 

AZ/DME 

L “ simple spreading loss is a function of the wavelength (X) and the 

slant range (SR): 

AZ/DME . , , 

^ PATH " 20*LOG^q (X/%SR^2/dmE^ 

K ppig - constant that includes transmitter power and circuit losses for 

the DME downlink (dBm) 

AZ GS 

G .j, >• azimuth ground antenna gain computed as a function of 


RW 

0 AZ/DME 


ORB AZ/DME 

G p - Orbiter antenna gain along Az/DME look angles (same 


ORB AZ/DME 

as G ^ and, hence, will be denoted by 

Q ORB AZ/DME) 

K AZ ” constant that includes transmitter power and circuit losses for 
the azimuth uplink 


EL GS 

G ,p - elevation ground antenna gain computed as a function of 


(dBi) 


(dBi) 


(dBm) 




RW 

EL 


ORB EL 
0 R 


Orbiter antenna gain along the Orbiter to elevation look 
angles 


(dBi) 


(dBi) 
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■ EL -AZ/DME 

L - same as L but using SRgj^ as the slant range 


(dBm) 


K el “ constant from elevation transmitter power and elevation uplink circuit 
losses (dBm) 


Equations for the above, and values to be used, are listed in section 10. 4. 3* 

RF lockon flags are then generated for each measurement (Az, El, and DME) by com- 
parison of the received signal levels to the assumed receiver thresholds in sec. 
10.4.3: 


If 

s 

DME 

DOWN - 

THRESHOLD 

DME 

If 

s 

AZ 

> 

THRESHOLD 

AZ 

If 

s 

EL 

> 

IHRESHOLD 

EL 


Then RF lock on is present -for DME, azimuth, and/or elevation, respectively. 

10 . 3.3 Data Validity 

It must be ascertained that both coverage constraints are met and RF lockon is 
present for each information channel before the data are flagged valid. 

10 . 3.4 Error Modeling 

Azimuth, elevation, and DME corrupted measurement formulations will be discussed 
separately. The dominant error sources are identified, briefly explained, and 
assigned to either the ground station, path, or navset. Those due to the ground 
station and path are assumed to be seen equally by all navsets. A particular 
hardware configuration has been assumed and is defined in the section on special 
notes . 
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Bias Md noise errors are identified separately. Due to the different types of 
biases present in eaoh measurement, different bias dependences are explained. 
Noise errors are assumed to be purely random with a normal distribution j how- 
ever, the variance may change as a function of known parameters. Different addi- 
tive bias formulations are suggested to account for the following types of bias 
mechanisms : 

a. Purely deterministic biases whose values would be the same the next time 
that trajectory is flown 

b. Calibration uncertainties idiose magnitudes may change from flight to flight 
but, having been randomly sampled, behave deterministically on that flight 
(ref. 10-3) 

c. Correlated biases whose exact values are unknown but whose effect may be 
modeled via an exponentially correlated random variable (ECRV) whose time 
constant may change (ref. 10-4), 

10 . 3 . 4.1 Azimuth Angle Errors 

The dominant bias sources are azimuth boresight error, angle pickoff bias, 
si delobe multipath, and beam truncation. Azimuth boresight error and angle 
pickoff bias are due to ground station alinement uncertainties. Therefore, 
a random number is chosen from a unit variance, zero mean, normal distribution 
to be multiplied by the formulation described in section 10.4.4. 

Multipath cannot be modeled deterministically due to lack of knowledge about 
reflecting surface, locations and sidelobes. The magnitude and rate of change of 
error is inferred from flight test data and modeled as an ECRV. 
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Beam truncation is purely deterministic and results from receiving only a por- 
tion of the scanning beam near the edge of the coverage volume. This effect 
begins at a true azimuth angle ('1^ TRUNC^ typically about +13.5° and continues 
out to the edge of the coverage volume (about +17°). An error is introduced 
equal to half the difference between AZ <rRUE TRUNC» always in the direction 

yielding a smaller corrupted azimuth angle. 

Of all the azimuth random noise sources considered (ref. 10-3)' beam pulse 
resolution error is the largest when signal levels are high, but at lower levels 
receiver noise bearashape error causes quite erratic measurements. A random num- 
ber is sampled from each measurement to be mulitiplied by the one-sigma random noist 
value computed as a function of signal level and geometry. 

All the bias may be assumed attributable to the ground station and path, but ran- 
dom noise is inserted on transmission and detection. 

The corrupted measurement formulation for azimuth becomes 

CAL TRUNC MP GND NAVSET 
^^CORROPTBD = *^TRUE ‘^BIAS '^BIAS '^BIAS %0ISE + '^NOISE 

where 

CAL 

T bias ” Calibration accuracy bias 
TRUNC 

"i bias Beam truncation bias 
MP 

^ BIAS - Multipath 
GRND 

NOISE “ Random noise due to the ground station 
NAVSET 

’J' noise - Ruudom noise due to the navset 
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10.3.4.2 Elevation Angle Errors 

Elevation angle bias sources are the same as for azimuth, but initial alinement 
uncertainties yield elevation axis level errors instead of azimuth boresight 
errors. Angle pickoff bias, multipath from sidelobes, and beam truncation all 
have similar contributions as with- azimuth. However, main lobe multipath occurs 
at true elevation angles below about 2 degrees, and will be assumed to behave in 
a purely deterministic fashion (ref. 10 - 5 ). 


Elevation angle random noise sources and dependences on signal level are about 
the same as with azimuth (slightly smaller in magnitude). It may again be 
assumed that biases would be seen the same in all navsets, but that random noise 
is inserted upon transmission and detection. 


CAL 

^^CORRUPTED - 4^TRUE + ®BIAS 


TRUNC 

®BIAS 


SLMP MLMP 
®BIAS + ®BIAS 


GRND NAVSET 
■*'®H0ISE + ®N0ISE 


where 


CAL 

9 bias “ Calibration accuracy bias 
TRUNC 

0 bias “ Beam truncation bias 
SLMP 

9 bias ” Sidelobe multipath 
MLMP 

0 BIAS " Mainlobe multipath 
GRND 

0 noise ~ Random noise due to the station 
NAVSET 

^ NOISE “ Random noise due to the navset 
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10. 3. ^.3 DME Errors 

As was the case with angle errors, three types of bias are present. The largest 
is due to receiver threshold mean times, idiich are expected offsets that are a 
function of signal strength. Incoming interrogation and reply pulses have a fi- 
nite rise time causing weaker signals to be detected later, resulting in a posi- 
tive error. The system is calibrated to minimize this error at the closer 
ranges (higher signal levels). A calibration bias terra is included tdiere a ran- 
dan number is multiplied times a one-sigma f emulation. The multipath may be 
modeled as an ECRV. 

Considering DME randcai noise, at high-signal levels the ground and airborne 
clock resolutions and airborne -output increment errors dominate with approxi- 
mately equal error contributions. At low levels the airborne and ground thermal 
noise become the largest. The formulation for the corrupted DME measurement is 
then: 

GRND OFFSET NAVSET OFFSET GRND CAL 
^^CORRUPTED = ®^TRUE ^^%IAS ^^BIAS ^^^BIAS 

NAVSET CAL MP GRND NAVSET 

+ DMEgj^3 + DMEgj^3 + DMEjjojgg + 0MEjjQj3g 

where 


GRND OFFSET 

BIAS “ expected bias 

NAVSET OFFSET 

DME Bj ^3 - Airborne expected offset 

GRND CAL 

DME 31^3 - Ground calibration uncertainty 

NAVSET CAL 

Bjj ^3 - Airborne calibration uncertainty 

MP 

DME - Multipath 
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GRND 

DME igoiSE “ Randan noise due to the station 
NAVSET 

DME jjQj 3 g - Random noise due to the navset 

10 . 3.5 Special Notes 

The algorithms presented and values to be used are subject to the following con- 
figuration and other assumptions 

a. DME double computes, viiich results in erroneous DME measurements and a 
required user software fix for ALT is assumed fixed in the hardware. 

b. Signal attenuation due to rain is not modeled. It could be easily inserted 
reflecting either rain cell or widespread rain conditions (ref. 10-1). Rain 
attenuation (dB) values are proportional to the propagation distance through 
the rain and the rainfall rate . 

c. Azimuth polarization errors are not modeled as yet due to lack of knowledge 
of the Orbiter antenna polarization characteristics. Although error magni- 
tudes are unknown, significant differences in the angle detected in differ- 
ent receivers are possible since three receiving antennas are used that may 
have differing polarization characteristics. 

d. Possible azimuth and- DME dropouts due to not having RF multipath fences are 
not modeled . 

e. Multipath error magnitudes reflect performance without RF fences present. 

f. A somewhat idealized Orbiter antenna pattern unaffected by the Orbiter nose 
boom present on OV-101 is assumed. 
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g. The effect of receiving only part of the scanning beam on signal levels at 
the edge of the angular coverage volume is neglected". ► 

h. Station location uncertainty, although traditionally modeled, is neglected 
due to expected survey accuracies. 
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10.4 MSBLS MODEL VALUES 

The contingent values for Orbiter and ground antenna locations, coverage 
constraints, RF parameter computations, and error model formulations are listec 
below. 


10.4.1 Geometric Observable Values (ref. 10-1) 

BODY 

- ORB ANT = Xq - 1.3 (ft) 

RW ^ ^ ' 

- AZ/DME = 16300 X py + 300 ypy - 7.6 Zp^ (ft) 

RW /s 

- EL “ ^RW — ^RW “ ^’6 Zrw (^6) 

The plus or minus sign depends on whether the installation is on the left or 
right of runway centerline. The minus 7.6 and 4.6 feet represents the height c 
the antenna above ground (plus Zp^ is down) . 

10.4.2 Coverage Check Values 


SR MAX = 119 090 

(ft) 

'f'MAX = 17-0 

(deg) 

^ ©MAX = 31.3 

(deg) 

4.3 RF Parameter Values 



10.4.3- 1 Ground Station Antenna Gains (ref. 10-2) 


DME GS RW 

p =12 + 10 sin (9(0 j^z/dME " 3-5)) 

RW RW 

- 20(l-cos(2.15*'t'AZ/DME>) ‘ -l^*® aZ/DME f^Bi) 
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If 1 '}' 

RW 

AZ/DME 

1 > 45® ; G 

DME GS 

R = -1° 

(dBi) 

AZ GS 
G ^ — 

5 

n=0 

RW 

AZ/DME^ 


(dBi) 

Where 

0 

0 

II 

2.10407996 

El 



Cl = 

4.84938878 

EO 



C2 = 

-8.06039523 

E-1 



C3 "" 

5.59289151 

E-2 

' 


C4 = 

-1.91706739 

E-3 



C5 = 

2.64114983 

E-5 



EL GS RW 

G ^ =27 cos (1.336*'}' EL^ 


(dBi) 


, RW , . EL GS 

If I '}' el I - 5 ^ T = 


(dBi) 


10.4.3-2 Orb iter Antenna Gains 


ORB AZ/DME BODY BODY . ORB EL 

G is a function of 0 j^z/DME AZ/DME ^ R 

BODY BODY ' 

be computed using 0 '}^ EL * Since the same antenna curve fit is used 

for both lines of sight the equation may be generalized to the form 

f „B0DY 2 ,„B0Dy2..., 

If (2.56 (0 - 2) + '}' p/^< 1|0 


ORB 

G 



n=0 


2.56 


„ BODY , BODY‘S 
(0 ■_ 2)2 + '}' 


]n/2 


+ 10 


(dBi) 
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Where 

ao = -5.578790329 

E-2 


= -1.916198777 

E-2 


a2 = 8.122955059 

E-4 


= -1.066190991 

E-4 
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ORB 

Otherwise G = -5 


(dBi) 


This low antenna gain represents exceeding the Orbiter antenna specified cover- 
age volume. 

10.4.3.3 Polarization Loss 

The polarization mismatch angle (4> pQ^) is assumed to be the angle between the 

Orbiter Z-axis and the runway X-Z plane. The orientation of the Orbiter Z-axis 
RW 

(z ORB^ is contained in the last column of the body to runway transformation 

matrix (T~^). <j> poL “lay be expressed as the complement 

RW EW 

of the angle between Z Qpp and the runway Y-axis ( Y ) . 

RW RW RW RW 

<i> POL = 90 - cos"^ (Zqrb’X /|ZorbI 1I (deg) 


The voltage interaction between the transmitted wave and receiving antenna, 
assuming linear polarization, is proportional to the cosine of (}> pQ^. The power 
loss in dB may be expressed as 

^ POL = 20*L0 G^q I cos (j) pQL [ (dB) 

10.4.3.4 Path Loss 

Assuming a wavelength (X) corresponding to a frequency of 15.7 GHz, the 
equations for path loss reduce to 
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^Z/DME 

L 

-46 - 20*L0G.|0 

(SR 

AZ/DME^ 

(dB) 

EL 

L = 

-46 - 20*LOG.|0 

(SR 

EL^ 

(dB) 

10.4.3.5 Constants 




The following 

transmit power 

and 

cable loss assumptions 

lead to four constants 

that are summed with the computed parameters above yielding the received powers 

Link 

Transmit 

power 

(dBm) 


Transmit 
cable loss 
(dB) 

Receive 
‘Cable loss 
(dB) 

DME down 

63 


- 3.8 

-1.7 

Azimuth up 

62 


-2.5 

- 3.8 

Elevation up 

62 


-2.5 

- 3.8 

^ DME DOWN 

= 57.5 


(dBm) 


K AZ 

= 55.7 


(dBm) 


K EL 

= 55.7 


(dBm) 



10 . 4 . 3.6 Receiver Thresholds (ref. 10.6) 


THRESHOLD 

DME = “77 

(dBm) 

THRESHOLD 

AZ = -7^ 

(dBm) 

THRESHOLD 

EL = -74 

(dBm), 


10.4.4 MSBLS Error Model Values 

The majority of the error formulations presented are adapted from vendor (AIL) 
theoretical error analyses (ref. 10-9 )• Exceptions are referenced. Flight test 
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data indicate that real world measurement histories behave as theoretically 
predicted in general. 

10.4.4.1 Azimuth Error Values 


w 3 r Ml ™ 

^ BIAS = Rni* 10“^* 202 + 113(1 + cos (2*^ AZ/DME) ) (deg) 

RW RW 

+ 225*sin2(2*'^ AZ/DME)*tan2(Q aZ/DME) 

¥ 

Where RN^^ = random number chosen from a unit variance, zero mean, normal dis- 
tribution at the start of each approach and landing to remain con- 
stant through rollout 


Y (Icr) varies between 0.0142° and 0.017°. Modeling simplifica- 

tion is achievable by assuming the conservative 0.017° 
(la) is constant (ref. 10-9). 


If AZ ^puE ^ 13*5° 


TRUNC 

^ BIAS ■ ■ (AZ^gyg-13.5)/2 (deg) 


If AZ true < -13.5° 


TRUNC 

^ BIAS "" - (AZtpue^‘' 3.5)/2 (deg) 


Otherwise T BIAS ~ 

The general formulation for an ECRV to model multipath for each measurement 
will be presented once. Thereafter, only correlation distance and one-sigma 
values to be used are listed (ref. 10-8). 
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First measurement: 


BIAS MP (t=0> - ^ MP * ^^2 


Thereafter: 


BIAS MP (L) = a*BIAS (t-^t) + cr f^p ( 1 - *RN 2 


MP 


Where Cf mp = sigma value 


(deg or ft 


RN 2 = randcm number chosen each measurement 


a = exp ( - At/Tgg) 


Tbs = DcoR *At/{SR(t)-SR(T~At)} 


SR = SR 


AZ/DME 


or SR 


EL 


D QQp = correlation distance 


(sec) 

(ft) 

(ft) 


For Y 


SLMP 

BIAS ■ 

Y 

SLMP = 

Y 

COR = 


.015 

6000 


(deg) 

(ft) 




GRND 

NOISE 




RW 


= RN2*10“3* ( 87 + 4.8*Y a 2/DME^ (deg) 




Where RN^ - random number chosen each measurement 


NAVSE-T 'J' , . , 

NOISE = RNi}*10"3» ^ 30 + 8.8*10 (-^-SAZ/10)}1/2 
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'H 

Where RN|^ - random number chosen each measurement 


10.4.4.2 Elevation Error Values 

CAL ® o r r 

®BIAS = 10-3{202+1113(1+cos(2%l)) 


+ l69sin2(20EL)]tan2('PEL)'}1/2 (deg) 


0 

Where RN - randan number chosen at the first of each flight and to 
1 Y 

remain constant (different from RN^^) ■ 

'CAL 

Note 0 BIAS assumed a constant 0.017 (icf) with little sacrifice 

V, 

in accuracy as with i utac (ref. 10-9). 


If EL true ^ 29 . 0 ° 

TRUNC 

® BIAS ■ “ (ELTBug-29.0)/2 (deg) 


TRUNC 

Otherwise 6 = 0 


(Note: Elevation beam truncation at the lower elevation angles is included in 

the main lobe multipath maximum value curve fits presented later.) 


SLMP 
® BIAS 

- 


0 



^ SLMP 

= .015 

(deg) 

0 



^ COR. 

= 6000 

(ft) 

MLMP 
® BIAS 

10-5) - 
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The path length traveled by the reflected ray (R ref^ roust be computed for 

RW 

for comparison with the direct ray path length (SR gj_,) in order to determine 

^LMP 

the phase difference between the two paths. © is then assumed to vary 

sinusoidally as a function of this phase angle between upper and lower limits 

RW 

which change as a function of elevation angle (EL XRUE^* ^ REF computed 

RW 

similarily to SR gg using +4.6 feet for the ground antenna location runway 
Z component. This represents the antenna image in the reflecting plane (the 
ground). The phase difference is then found by taking the fractional portion 
of the path difference divided by the number of feet per wavelength and multiplied 
by 360 ®: 


^ MLML ^ RW RW 

©el = 360*FRAC {(Rref - SEgg)/ .06382} (deg) 


Maximum and minimum error bounds are described by third order curve fits: 


© 


MAX 

MLMP 


3 

S 

n=0 


^n * 


RW 

EL 


n 


0 OT ) 


(deg) 


MIN 3 n 

® MLMP = 2 - bjj * ( 0 el ) 

n=0 


fl/here 


Sq = 0.3786178992 

bo = -0.5041173814 

a-j = -.-5765868845 

b.] = 1.599867284 

a2 = .3409840709 

b2 = - 1.231612944 

33 = -. 0700743883 . 

63 = .2834967565 


(deg) 
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The size of the mainlobe multipath error envelope (A,) is then, 

MAX _ MIN 
■^ = ® MLMP ■ ® MLMP 


With a mean value (B) of 
MAX MIN 

® ^ ® MLMP ® MLMP^''^ 


MLMP 

The previously computed phase angle is used to determine 0 


If EL .jigug < 2° 

„ MLMP . MLMP 

® BIAS = * (0 2L ^ ^ 

2 


MLMP 

Otherwise 0 = 0 


GRND 0 ^ , RW 

® NOISE = RN3*10“3 . { 30 + 1.6(0 gg) (deg) 

NAVSET ' 0 

® NOISE ~ RNi|*10~^* ■f 30 + 3*3*10^“^’'^EL''^^®^}^'^^ (deg) 

10.4.4.3 DME Error Values 
GRND OFFSET 

DME = 0. 0068*10" ^^DMEDOWN/20) .,3.7. (ft) 

NAVSET OFFSET 

OME gjj^g = 18.0 (ft) 


GRND CAL DME 

DME 3 j ^3 = ^ I * 5.0 (ft) 
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DME 


NAVSET CAL 
BIAS 


RW ^ 

= RN 2 "" * { 71 + (1.'97 + 0. 00005 *SR^2/DME^ ^ 


DME 


1/2 


For DME 
DME 

^ MP = 
^ COR = 


MP 

BIAS “ 
10 

6000 


(ft) 

(ft) 


GRND DME /aq 1 /^ 

DME noise = RN 4 » { 25 + 0.15*10 ^"°"^DME DOWN-^^°^ } (ft) 


NAVSET DME ^ „ /1n^ i 1/? 

DME noise = RN 5 * t 27 + 0. 15*10 *‘“°"^DME DOWN^^'J^ > (ft) 


DME DME DME 

Note: RN and RN 2 3 .re chosen once per flight while RN ^ (for 


DME IWE 

ECRV) , RN 4 , and RN g are sampled each DME measurement. 
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APPENDIX A 
REVIEW OF CONCEPTS 

A. 1 RANDOM PROCESSES 

A random, or stochastic process, is a set of random numbers each representing, 
in some way, the probabilistic outcome of a particular event. It is not neces- 
sary, in' the understanding of the error models contained in this document, that 
one mathematically perceives the meaning of random processes since only a few 
simple types are used. The most simple, random noise and fixed bias, involve the 
multiplication of a representative standard deviation and a random number chosen 
from a given distribution. This section will discuss the methods by which a dig- 
ital computer is utilized in generating random samples from two such distribu- 
tions: (1) uniform over the interval zero to one, and (2) Gaussian (normal) 

with zero mean and unit variance. 

A . 1 . 1 PSEUDORANDOM NUMBERS 

Chambers (ref, A-1) presents an excellent- treatise on the generation of random num- 
bers on digital computers. The reader is referred to this paper for more com- 
plete detail of the following description of the particular technique being pres- 
ently used for error analysis studies. 

The general form of the random number generating function used on digital 
machines is 

= ax^ + c (mod m) 

4 * 

where xj = j random number 

a = a multiplicative factor 

. ^ ■ A-3 
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c = an additive factor 
m = a positive integer larger than Xj 

This form is called the mixed congruential method. The special case where c = 

0 is called the multiplicative congruential method. The behavior of the tech- 
nique in terms of its period (quantity of random numbers before repetition) and 
its randomness (uniformity of distribution of consecutive pairs) is determined 
by the choice of a and c. It has been shown that a multiplicative congruentia] 
method (c = 0) with the multiplier, a, on the order of Vm' provides excellent 
randomness. Note, also, that the recursion relationship is deterministic. 

That is, the sequence of random numbers is repeatable for a given starting 
number, x^; hence the term pseudorandom numbers. 

The technique used on the Uni vac 1108, as programed in subroutine RDM, is a 
multiplicative method with a = 5^^. The choice of this constant has been based 
on randomness, it would seem, since a number based on speed of execution would 
be an integer power of 2 (allowing the multiplication to be reduced to bit 
shifting). The modulus, m, is chosen as the largest positive integer possible 
in the 36-bit word (377, 777, 77-7, 777s). practice, this is merely the least 

significant 35 bits of the word in question. Actually, the multiplication 
by a causes (eventually) the integer to overflow into the adjacent register. 
When this occurs, the less significant register (that containing the less signif- 
icant numbers) is used to extract the random number. Once this occurs, the inte- 
ger bit pattern is converted to a floating point number and the mantissa (least 
significant 27 bits) is taken to be a uniformly-distributed random number on the 
unit interval. 
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Mathematically, it looks like this. 

Xi+1 = (mod 377,777,777,7778) 

Ui +1 = x*+i (mod 1) 

where u is a uniformly-distributed random number on the unit interval. 

A. 1.2 OTHER DISTRIBUTIONS 
Relatively few random processes are uniformly distributed, and it therefore be- 
comes necessary to convert a uniform distribution to one of another distribu- 
tion. The procedure is as follows. 

let Py(u) = 1 0 < u < 1 (uniform) 

where p is the probability density function. 

It is desired that y = f(u) has a density of p(y). 

If we can calculate 

U = f’'(y) 

then 

dP = Py(u)du = 1 du = p(y)dy 
Therefore , 

u = f”^(y) = ^(y)dy + c 

(where c is a constant of integration) is the desired transformation. 

( 

*Note that the integer has been correlated to a floating point number. 
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A. 1.3 GAUSSIAN DISTRIBUTION 

The Gaussian, or normal, distribution is most frequently used in error analysis 
and is generated, through an intermediate process, in the following manner. 

let Pu(u) =1 0 < u ^ uniform 

Produce a Rayleigh-distributed variable r by 

r = a(-2 In u)^/2 

where the probability density function for .this variable is 


p(r) 



r 


2 


20 


2 


e 


It is now possible to produce a pair of correlated Gaussian variables using 
r and a tmiformly-distributed variable 9: 

X.J = r COS0 
X 2 = r sin0 


where 0 is uniform over zero to 2 it. 

Substituting and constraining the distribution to have unit variance, 

XI = (-2 In u-j )^‘^^cos27TU2 
X 2 = (-2 In u-|) ^^^sin 2 fru 2 

Proof of this method in its ability to produce Gaussian variables can be 
generated by an inverse procedure. 


A-6 



ORIGH\?AL PAGE IS 

OF POOR QUALETY 79FM5 

Assume and X2 are independent Gaussian variables with densities p(x^ ) 

and p(x2)* Then, the joint density is 


p(x-| , X2> = p(x-] )p(x2) 
1 


,- 1/2 


(X-]^ + 


2 iTa^ 


r 2 

1 20^ 

p{x^, X2) = e 

2 ira^ 

pop 

where r = x-j + X2 

To complete the substitution of variables, assume the new probability density 
to be of independent variables r and 6 . That is 

_ ^ 

■] 2 ct^ 

fCj^j e) = e 

2 iTa^ 

By differential probability 

dP(r, e) = f(r, e)r dr d© 
where P = cumulative distribution function 


0 < r < 00 

0 < 0 < 2ir 
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Rayleigh distribution 


1 

P(0) = — uniform over 0 •* 2-n. 

2tt 

This process is programed in subroutine RN2S. 

A. 2 GAUSS-MARKOV PROCESSES 

A Markov process may best be explained by conditional probabilities through a 
) 

I 

simple exampl^ (ref. A-2). Let Ep E^, . . . . , E j be a complete system of events 
Consider the outcome of a sequence of consecutive statistical trials (flipping 
a coin, rolling a die, etc.) from the view of the occurrence of the events 
Ej (j z 1, 2 ...) and define the random variable ^ follows: 

Cn = j if Ej is the outcome of the n^^ trial. 

When all trials are statistically independent 

P{?n = jiCo = io> = il ••• Cn-1 = Vll = P(?n = 

That is, the probability that the event will occur on the n^^ trial is 
independent of the outcome of all previous trials. 
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If we define a sequence where the outcome of the trial depends only on 
the preceding outcome (n-1), we have defined a first-order Markov process. 

That is, 


P{?n " jUo “ ^o’ ••• ?n-1 “ ^n-l} " ^I?n " ^l^n-l ~ ^n-ll 

I 

I 

Extension of this definition to dependence on p prior realizations defines 

th * 

a p -order Markov process. The remainder of this discussion will concern 

itself with processes of order one. 

It is reasonable to assume 


E(w(k)) = 0 


E(w(k) w(j)) = 


lc=j 

k;^j 


E(X(o)) = X(o) 
E((X(o) - = Po 


E(X(o) - pjj) w(k)] = 0 for all k 


(pure noise) 


(unbiased) 


If the statistics of' the random variables are Gaussian, the process is called 
^ Gauss-Markov. A scalar Gauss-Markov sequence and the variance of that sequence 
may be represented as (ref. A-3) - 


X(k+1) = a(k+1, k) X(k) + w(k) 
P(k+1) = a^(k+1, k) P(k) + q 
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where 


X(i) .= random variable at i 
ct = correlation coefficient. 
w(k) = purely random noise at k 
P{j) = variance of X(j) 
q = variance of w 


QF 


I 

A solution for X(k+1) and P{k+1) in terms of X(o) and P(o) may be found 
through their recursion relationships by 


X(k+1) = X(o) + w(o) + w(1) + ... + w(k) 


or 


X(k+1) = X(o) + ^ w(j) 

j=0 


and 


P(k+1) = P(o) + q(o) + q(1) + ... q(k) 


but for 


q constant 
P(k+1) = 


k 

P(o) + q 2 

j=1 


If a^<i , 

P(k+1) = p(o) + q 

For the process to be stationary, P will approach a limit for k-*-«>, which 
can be satisfied for a^<1. Taking the limit, 
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lim, P(k+1 ) 
kx» 


lim 

k-x» 


o: 


2k 


P(o) + q 


1-a' 


2k 
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and 

q 

P(k) -> — 
1-a^ 


The standard deviation of w(k) = Vq" allowing w(k) = u(k) where u(k) 
is a unit-variance, zero-mean Gaussian random variable 


Now, 


X(k+1) = a X(k) + V 



or 

X(k+1) = a X(k) + u(k) 

where = standard deviation of X(k) , 

which is. the desired form. 


For exponentially-correlated, time- related variables 


a 


e-AT/T 


AT = the time between k and k+1 


T = correlation time constant 


causing the -form of the process to become 

X(tjj+-)) = X(tij) + -^i-e“2AT/T 
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A. 3 QUANTIZATION ERRORS 

Quantization, as the name implies, connotes the value of parameter being 
truncated or otherwise rounded off to some level of significance. A simple exam- 
ple of this concept is the digital computer where word length constrains the in- 
ternal numerical significance and transformations between internal and I/O num- 
ber bases cause further numerical inaccuracies . 

For navigation systems the quantization of data due to the physical devices 
(both hardware and software) may have one of two effects. The first, which we 
might call destructive quantization, truncates or rounds the data to a certain 
level of significance and forever loses the knowledge of the remainder. Mathe- 
matically this appears as 

^actual “ ^actual “ ^ actual 9 

where q is the quantization level 

A second type of quantization occurs- when the actual system is retaining the com- 
plete significance of the data but the device utilized to read or sense the data 
is passive (does not disturb the data) and has a higher quantization level than 
the raw data. A simple example of this is a time clock used for punching times 
on cards, letters, etc. Although the time appears in truncated form (nearest 
minute, nearest second, or whatever)' the internal clock is keeping time to a 
much higher degree. That is, the data not shown in a given output will be 
retained and manifested in the next output. 

An example of this nondestructive quantization in navigation systems is a pulse 
accelerometer. This device accumulates sensed delta-velocity until a particular’ 
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threshold is reached where the accelerometer sends a pulse representing a partic- 
ular quantum of output. So at .any one time the observer has only the largest 
integer pulses available to him, but the device has retained the remainder. 

This remainder, of course, is used as the initial accumulation for the next 
pulse. Mathematically, 

^sensed ~ ^actual ~ ^actual 
Note that the actual value has remained imperturbed. 

Nondestruct quantization results in a correlation of the error in a given 
quantity. That is, the error or uncertainty in X(i) is dependent on the 
error in X(i-1). To return to our time clock example, if A punches in at 
8:00:59 on a clock that truncates to the last minute, his card will read 
8:00 and will be 59 seconds in error. B punches in 2 seconds later at 8:01:01. 
His card reads 8:01, which is only 1 second in error, but he gets his pay 
docked anyway. Even though they arrived 2 seconds apart, they appear to 
have arrived 1 minute apart. 

Let's look at the pulse accelerometer again to determine the error for destruct 
and nondestruct quantization. Lear (ref. A-4) presents the following example. 
Assuming a gravity-free environment and linear motion, consider nondestruct 
data. 


A 


A 


^i = %-1 


6Vi 


where 6V^ = + q^ - - q^_., 

A 

At the end of the first sample period the velocity estimate, R^, is 


A-13 



OF POOR QUAUTY 


79FM5 


Rl - Rq + q-| - Vq - q^) 

= + 9l - <3o 

for R^ = true velocity at time 1. 

At the end of the second period , 

A 

R2 = R2 + Q 2 - Qo 
^ v» 

and expanding to the i^“ period, 

A 

Ri = Ri + qi - q^ 

The standard deviation, then, of the velocity error due to nondestruct quant- 
ization is 

Og = constant 
R 

For the case of destruct accelerometer data the change in velocity over the 
interval i-1 ^ i is 

6Vi = Vi - Vi_i + ,i 

Note that the error in is uncorrelated (dependent only on the quantiza- 

tion error at i, q^^). Performing a similar analysis as for the nondestruct 
case 

A 

Rj^ = R^ + (qj^ + q^_-] + ... + q^) 

Calculating the standard deviation. 

Or = yr Oq 
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where i is the total number of samples. 

Note that the error statistics are not constant and represent a larger error 
than nondestructive data for i>2. 


A. 4 COVARIANCE MATRIX 

Most error analyses, linear or Monte Carlo, utilize covariance' matrices in sev- 
eral ways. Those most important to our discussion are the input and output 
matrices used for initializing and 'displaying results, respectively. The follow- 
ing discussion will assume a basic knowledge of statistics with respect to 
multivariate distributions, statistical moments, etc. 

Given the n-dimenslonal vector random variable x with joint density function 
f(x) we can define the mean and variance of each component as well as joint 
or mixed moments in the following manner. 

A. 4.1 MEAN OF i^^ COMPONENT 

CO 

Ui = e(x^) = J' x^fj_ (x) dX£ 


where f^ = marginal density function of x^ 
E = expected value operator 

A, 4. 2 VARIANCE OF i^^ COMPONENT 


= E[(Xi - ^i)^] 


(xi - ^’i(x) dxi 
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Oij = e((xi - ViKxj - Pj)) 



(x^ - iii)(xj - tij.) (x)dXj^ dXj 


The mixed moment, a,- j is sometimes called the covariant term and is actually 
composed of the following elements. 


^±3 - Pij ^i ^3 


where = correlation coefficient (linear) 

Ojj = standard deviation of the variable, xj^ 


A. 4. 4 CORRELATION COEFFICIENT 

The variable p,- j above is shown to obtain the values 

-L J 




1 


by proving the following inequality 



< O-; a. 


utilizing a technique outlined in reference A-5. 

Consider the real variable X and the function of X 
f(X) = E((xi ~ Ui) - X(xj - Uj)}^ > 0. 

It must be true that f(X) > 0 since it represents the variance (by definition) 
of a zero-mean, random variable for all X. Expanding 

f(X) = 0^2 _ 20ij X + X2 > 0 
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It follows that the quadratic f(^-) has either no real roots or two equal real 
roots. The conditions for this is that the discriminant satisfies the inequality 

“13^ - “ 3 ^ ; 0 

or 

CT . -2 < O ' . 2 0-2 

If there are no real roots, f(X) > 0 and we have 
°13^ < a/ 

Suppose there is a real number of such that f(X^j) = 0, it would follow 

that that value of we would have 

z = Cx^ - Pi) - - IJj) = 0 

for f(Xg) = e(z^) = 0 

since the variance of a random variable is zero if (and only if) that variable 
is equal to its mean with probability one. Therefore, for the following equality 




2 



a 


.2 

J 


holds only if there exists a linear functional relationship between Xjj^ and 
X4 as defined as z above. We define this relationship to be 


“13 = “3 


such that 


-1 < p . . < 
- - 


c-l> 
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The linear ftinctional relation is, of course, 
poi 

^i “ ^i “ p = +1 . 

We may loosely presume the correlation coefficient , p j , to be some measure of 
the tendency of the two variables to be linearly related (that is, the tendency 
of samples of the variables to cluster about a line). This would imply, in esti- 
mation theory, the observability of certain parameters; that is, the 
separability of the effects of each within the observations. This knowledge 
could be very useful in predicting the outcome of a given experiment. We will 
use, as an example, a simple least-squares problem. 

Given that we are trying to predict a value for x, namely x, that minimizes 
the sum of the squared deviations, min e( (x - x)^), we note that this yields 

/N 

X = Suppose we have joint random variables (x, y) and we wish to again 

predict x by observing y (linearly) . 

We therefore seek 

X = ay + b 

vftiich will minimize 

e((x - x)2j = eC(x - ay - b)^] 

This equation may be written as 

e[(x - Ux) - a(y - Uy) - b + - aUy]^ 
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and is minimized by 


. b = - ayy 

CTxy PCJx 



The best estimate becomes 


X 


P^x 


(y - yy) 


Px 


Note that for zero correlation, p = 0, the best estimate is the mean, 

X = y^. 

A. 4. 5 COVARIANCE MATRICES 

Alluding again to our vector random variable x = (x-| , X2, ... Xj^) with joint 
density function f(x) = f (x^ , X2, x^), we may recall the following 

quantities 

yj^ = e(x^] = mean of 

= E((xi - y^)^) = variance of Xj_ 

^ij = Pij ^i ^j = e((x^ - yi)(xj - yj)) = 

covariance of x^^ and Xj 
= correlation coefficient between x^ and Xj 

It follows that 
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Another convenient form of the correlation matrix is 



The relationship of the statistics, and c?y, through use of the eigen- 

vectors and eigenvalues, is developed in figure A-1 geometrically and for two 
dimensions. The extension to higher dimensions should be trivial since we only 
consider two variables at a time. 


The covariance matrix shown above (A) as an equiprobability ellipse can be 
represented as the covariance matrix Ag in its eigenspace (^i, ^ 2 ^ t which 
case it would be diagonal, and as the matrix A^ in the space (x-|, ^ 2 ) ‘ The 
characteristic equation for the surface may be computed as (ref. A-6). 
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= x^ A"'' X 

where X is any arbitrary vector. 

If we chose to represent the matrix as A^ chose an eigenvector of 

length £2 X, then the above equation reduces to 

k^= 1 

Since this must hold true over all orthogonal transforms, 
k^ = Ax“^ z 

then 

2 ~ 

k = 1 for all z that intersect the ellipse 
and we will choose the vector a, o in x-space, which represents the intersect 

9 

of the ellipse and the x^-axis. The above calculation results in 
2 

k2 = a2/[o^^ (1,- p2)J = 1 

or 


' a = ox^ (1 - 

where it can be shown that 0 ^^^ is the intersect of a perpendicular to the ex- 
treme point of the ellipse with the x-|-axis (fig. A-1). Inspection shows 
that P 0 JJ completes the right triangle of side a and hypotenuse q^j.. The di- 
rection of that side depends on 0 and q ^ (the eigenvalues of a one-sigma 

el e 2 

ellipse) and 0 , the angle between g-| and x-]. 
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If we consider the matrix and rotate to x-space, and get 


” 



2 




cosf* 

sin© 


Op ‘ 0 


cos© -sin© 

c 

-sin 6 



2 



= 

cos 6 


0 Op2 


sin© cos© 


_ 

_ 


_ 


_ _ 

L 


(v"4) 




sin^P + Op cos^P 
2 


The off-diagonal term may be equated to its x-space notation 


I sin 26 (o^^ - Op^) = p^y O^ Oy 


resulting in 


Pxy <"x 


sin 


The parameter Oy has the same geometric interpretation as the extreme 

value of the ellipse. 


A . 4 . 6 ADDITIONAL PROPERTIES 

The covariance matrix is symmetric by its definition. It is also positive 
semi definite , as is shown by the following proof. 

Consider the vector constant a and construct the scalar random variable 
z = a^ (X - 

Then 


e[z] = E(a'^ (x - e[x - Px^ = 0 


and 
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= e(z^) = E[a'^(x - - Px)"^) 

= a^ e((x - Px)(x - Px)'*’)a 
= Ax a > 0 for all a 

since variances are always nonnegative. Aj^. is positive semi definite. For 
z=0 there exists a nonzero vector a^ such that 
_T - _ 

Zq = (x - Pjj) = 0 



indicating that there exists a linear functional relationship between the 
components of x (analogous to [p] =1). There are (n - k) such linear 
relationships where k is the rank of A^. 

Diamant (ref. A-7) has produced a parameter that is indicative of the elongation 
or eccentricity, of an n-dimensional ellipsoid (representing an n-dimensional 
covariance matrix) , which is called asphericity, and is computed as 
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The correlated random vector is computed from the mean and covariance matrix 
by 

X = vi + C a 

where C = a lower triangular matrix such that cc'^ = A 

U = a vector of zero-mean, unit-variance random numbers generated 
by the technique described in section A.1. 

This method is shown to provide the correct statistics by generating the 
covariance matrix for X and comparing that to A by (ref. A-8). 

(x - U)(x - = C U 

Taking the expected value of both sides 

e[(x - u)(x - u)^] = c e[u iF) c'^ 

but since 

e(.u u"^) = I 

Ajj = e[(x - y)(x - y)^} = C = A 

The matrix C is determined from A by the square root method (ref. A-9) 
as follows. 
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°11 


c 


/ 

i-1 

°ii 

= Uii - 

I 


\ 

k=1 


/ 

3-1 


= 

I 

k=1 

°ij 

= 0 i<j 



1/2 


'ik , 


i>1 


•ik ° jk I jj 


The distribution of the output vector will, of course, depend on the distribu- 
tion of the random vector U. Computation of the elements c. ■ is carried 
out successively by columns. 


A. 4. 8 OUTPUT COVARIANCE MATRICES 

Traditional output of error analyses is the mean and covariance matrix of 
the output vector random variable y - ]iy and Ay, respectively.' An alterna- 
tive method of interpreting error analyses results is presented below. 

A. 4. 9 DIAGONAL ELEMENTS 

Since we are all now experts on the properties of a covariance matrix, we 
immediately are drawn to the diagonal elements that we know represent variances 
of the individual components of the random vector. We are aware, however, 
that they represent marginal variances. That is, the variance of the i^^ 
component of x is; 




—CO 



Cxi - f(xi, X2, 


. .., Xjj) dx-j dX 2 ... dXjj 
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It may be interpreted as the variance in the i^h component given the statistics 
of the other n-1 components and the dynamics of the system. 

We may perform a similar analysis on several components at once when we realize 
that most of the vector variables we deal with contain components of different 
units (position, velocity, angles, etc.). The adding of variances with similar 
units produces a variance of the vector uncertainty magnitude- in a multi- 
dimensional quantity. That is, the uncertainty in the magnitude of the error, 
not the uncertainty in the magnitude of the quantity. Since elements with like 
units are generally grouped together, the added variances constitute a serai- 
trace of the matrix. The square root of this quantity is called the root semi- 
trace (RST) and represents the standard deviation of the quantity. Note that 
this produces the same number as the root sum square (RSS) of the standard devia- 
tions of the same parameters (a very popular quantity). 

A.il.10 OFF-DIAGONAL ELEMENTS 

Recall that an off-diagonal element of a covariance matrix is composed of 
a linear correlation coefficient and two standard deviations 

^ij = Pij ^i 

Correlation coefficients, as we have noted above, indicate the tendency toward 
linear relationships between the variables. Even though it may not appear so at 
first glance, a correlation coefficient of one indicates a functional linear 
relationship between two variables independent of their units. This implies 
that such a coefficient appearing anywhere in the matrix denotes a singular 
matrix. A simple proof of this is offered by W. W. Ferry. Given that 
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then the ratio of terms in the i^^ and rows is 

'^ik Pik Cfi cfk Pik Cfi 

cfjk Pjk Cfk Pjk cfj 

Now, it follows that if = +1 , then 

Xj^ - = c{Xj - ]ij) (linear dependence) 

and G is found by 

E[(x^ - - Uj)] 

Pij = = 1 

Oi Oj 

Oi Oj 0 ± Qj 

c^i 

c = — 

Again, constructing the ratio terms, 

*^ik Pik °'i ° s((Xj - yj)(xjj - pj^)) ffjL 

^jk Pjk E[(^j “ yj)(Xk - Pk)] 0j 

f^ik 

= — = constant 

« « 2 
^jk CTj 

The above proof implies that now j is a constant times row i 
condition for singularity. 


a sufficient 
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The implication of a singular matrix depends on the problem being analyzed. In 
estimation problems singularity suggests nonobservability or ill-conditioning 
due to the dynamics of the system. 

A. 4. 11 EIGENVECTORS AND EIGENVALUES 

Since a nonsingular covariance matrix is symmetric and positive definite, it 
will have a set of n-positive eigenvalues and corresponding eigenvectors 

U^. The eigenvectors describe the direction of the principal axes of an 
n-dimensional ellipsoid and are an orthonormal set. The corresponding eigen- 
values are the dimensions of these axes and may be used to represent the 
standard deviations of a set of independent random variables whose statistics 
are depicted in the original matrix (rotated into the x-space). That is, 



Insight into many problems may be gained by inspection of eigenvalues. A 
discussion of this, subject is necessarily waived except for a simple example 
(ref. A-10). Consider the covariance matrix P associated with the particular 
estimation process . 

y = Ax + rt 

where y is a vector of independent observables 
X is the true state 
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T] is a vector of noise 
A is a linear matrix operator 

^ 

where the best estimate of x, x may be 
X = {A^WA)-"’ A*^ Wy 
for W = a weighting matrix 

The covariance matrix associated with the error in the estimate is 
P = (a"^ WA)"’’ 

T 

Suppose that we cannot invert A WA for some reason (a criterion for deciding 
this ahead of time will be discussed later) and we wish to know what is wrong. 

We generate the eigenvectors u^ and their corresponding eigenvalues 
which we know satisfy the equation 

U = UA 

where U = a matrix of eigenvectors [u^ U 2 
A = diagonal matrix of eigenvalues 

We may diagonalize P“^ by 

UP""* = D 

where the components of D, d^j_ represent statistics of n-independent random 
variables and represent the information in this variable gained by observing 
the y- vector . 
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Suppose djj is small. This means that little Information can be gained 
about this parameter from our choice of observables and it would be best 
to ignore it or consider it another way. 

On the other hand, say dj^ is large, implying good determination of this 

variable. Now, if u^ (the eigenvector associated with the k parameter) 
had two large components (about equal) and the rest negligible, it is possible 
to determine a linear combination of the large components even if they are 
nonseparable. Insight into many similar problems may be gained by inspection 
of the eigenvalues and vectors of an information or covariance matrix. 

Pace has pointed out that one criterion for inversion of a matrix on a digital 
computer is that the exponent of 10 in the ratio of the determinant to the 
trace have an absolute value less than the number of significant decimal 
digits contained in the machine word. That is, for 

det(j:] 

= 10 “ 

tr(e] 

jm| < number of decimal digits 
is a necessary condition for invertability • 

A. 4. 12 ONE-SIGMA VECTORS 

By this time we should be familiar enough with scalar random variables to 
realize that their standard deviations represent a measure of the probability 
that a sample of one-sigma or less will occur. Similarly, we can assess 
probabilities to this many or that many sigma events. There does not appear, 
however, to be a simple extension to multidimensional random variables (random 
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vectors) of this concept. Detchmendy (ref. A-6) with reference to Cramer (ref. 
A-11), treats this subject excellently in a manner that is summarized below. 

Let X be a zero-mean, Gaussian random vector with covariance matrix A. 
Surfaces of equal probability are defined by 

xT A-1 ^ = k2 

where k = constant. This may be interpreted as any vector sample taken from A 
that satisfies a particular value of is equally probable. 


Now the probability that a sample vector of dimension m lies within the 
ellipsoid defined by k is 


m 

PCk) = 

m 

2111/2 r{- + 1) 
2 



j^m-l g-1/2 r2 


where r(i?.) = 



gamma function 


For the special cases of m = 2, 3> and 6, 
m=2 

P(k) = 1 - e-k^/2 

m=3 k 

2 -1/2j,2 

P(k) = - -ke + e dr 
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P(k) = 1 - _ 
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A plot of these functions is given in figure A-2. 

We note two interesting points: (1) any eigenvector u^^ of length is 

a one-simga vector (k=1), and (2) the length of a vector may not be repre- 
sentative of its probability of occurrence. That is, a short vector in the 
direction of the smallest principal axis may be less likely than a relatively 
long vector in the direction of the longest axis. 



Even though b is longer than a, it is more likely to occur. 

A. 4. 13 ACCURACIES OF SAMPLE COVARIANCE MATRICES 

A question that necessarily arises when generating covariance matrices by 
Monte Carlo Analysis concerns the accuracy with which the sample variances and 
correlation coefficients have been constructed in N samples. In the discussion 
below we will note two features: (1) the accuracy in constructing variances is 

independent of the size of the variance whereas, (2) the accuracy in con- 
structing a correlation coefficient depends on its sample value: 

p 

The accuracy with which the sample value of a variance S is computed in N 
samples (when compared to' its theoretical value c^) is given in figure A-3 (ref. 
A- 12) for various levels of confidence. At 90 percent confidence we see that for 
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30 samples the true standard deviation will be within +0.23a of the sample- 
value S. If we up the sample size to 100, this difference becomes +0. 12; 
not quite halving the error for more than tripling the sample size. Note 
the following. 

a. The value — is independent of S. 

a2 

b. Due to the flattening of the curves for increasing N, one must increase the 
sample size many- fold for lesser increases in accuracy. 

c. Extension of these curves will show that 1 percent accuracy in a can be 
achieved at 90 percent confidence for N around 550. 

d. The sample size most appropriate for a given problem depends on what interpre- 
tation is placed on the answers. If one wishes to assign probabilities to 
these standard deviations, he must know the distribution, of the results. 
Traditionally, we assume the results to be Gaussian (by the Central Limit 
Theorem)' when they are only approaching a Gaussian distribution. Such an as- 
sumption may completely invalidate the use of large sample- sizes. 

Correlation coefficients behave a little differently (ref. A-13). Define the 
foil owing terras . 

• 

r = computed value of correlation coefficient 
p = theoretical value of correlation coefficient 
N = number samples to compute r 

1 . I+P 

= - In 

2 1-p 
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Z has a nearly Gaussian 


1 1+p 

mean = - In 

2 1-p 
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If the distribution of the samples is Gaussian, then 
probability density function with 


1 

variance = — 

N-3 

The confidence limits for r may be determined by the following examples. 

a. Example 1.- Suppose for N = 100 observations a value of r = 0.8 is 
computed, then Z = 1.098 and the standard deviation is 0.101. Since Z is 
nearly Gaussian, then it will be within +1.65 a of its mean 90 percent of 
the time (i.e., the 90 percent confidence level). Thus, with 90 percent 

1 +p 

confidence 1/2 In — is within 1.098, +1.65 (0.101); i.e., within 

1-p 

0.931 to 1.265. Translating this to the correlation coefficient, p, then, 
with 90 percent confidence is within p^ = 0.?0 to p 2 = 0.852. 

b. Example 2.- In order to illustrate the manner in which the confidence level 
values for p vary with the computed value of r, consider the identical 
case as above except with r = 0.3. For a 90 percent confidence level 

on p, the interval works out to be p.| = 0.143 to p 2 = 0.40. 

So we see that as the correlation coefficient estimate r gets smaller its 
accuracy becomes less. But, as the coefficient becomes smaller, the results 
are less influenced by its size and the two effects cancel each other. 
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APPENDIX B.1 

COORDINATE TRANSFORMATIONS 

The coordinate transformation of primary interest is that between the reference 
inertial coordinate system and the (rotating) Earth-fixed geodetic coordinate 
system. For Space Shuttle work, the mean of 1950.0 coordinate system has been 
adopted as the reference inertial system. 

Numerous subroutines that exactly perform this transformation based upon the 
adopted values of the relevant astronomical constants are available. These 
subroutines are consistent to 14 or more significant figures with the adopted 
constants. Two such subroutines are documented in the following reports: the 

Houston Operations Predictor /Estimator (HOPE) Engineering Manual (ref. B.1-1) 
and Description of a Self-Contained Subroutine which Analytically Generates 
Interplanetary Coordinate System Transformations Referenced to the Mean of 
1950.0 Epoch (ref. B.1-2). 

The exact subroutines are much too long for use in an onboard computer. A 
simpler subroutine more suitable for onboard use is documented in volume II of 
the Space Shuttle Navigation Software (ref. B.1-3) documentation. 

This simpler subroutine is consistent with the more precise subroutines to about 
7-8 significant digits, or about 1 foot at the surface of the Earth. 

Errors in the transformation from the inertial system to the Earth-fixed system 
are generally ignored as a separate error source. 

A second coordinate transformation of interest is that between the Earth-fixed 
coordinate system and the topocentric coordinate system at a point on the Earth's 
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surface. The topocentric system is oriented with its axes in the vertical, 
east, and north directions. Equations defining the transformation may be found 
in reference B.1-3* 

Errors in this transformation are negligible compared to errors in quantities 
measured in this coordinate system (e.g., elevation angle). 

The various coordinate systems that will be used during the Space Shuttle pro- 
gram are described in the Coordinate Systems for the Space Shuttle document 
(ref. B.1-4). 
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APPENDIX B.2 
GRAVITY MODEL 


B.2.1 GENERAL DESCRIPTION 

The motion due to the gravitational interaction of two point masses is 
completely defined by the inverse square law of Newton. Newton's law also 
applies to radially homogeneous spheres. However, the gravitational potential 
of most physical masses found in the universe must be altered somewhat from the 
idealized two-point representation due to nonsphericity, nonhomogenity, 
rotational dynamics, and nonrigid deforraational effects. Methods of modeling 
the Earth's geopotential include spherical harmonic expansions, point mass 
representations, ellipsoidal harmonic expansions, density layer models, etc. 

An Earth gravity model formation due to Pines (ref. B.2-1 and B.2-2) requires 
that the orientation of a point in space be expressed in terms of direction 
cosines s = x/r, t = y/r, and u = z/r, where x, y, z are the Earth-fixed 
Cartesian components of the position vector of magnitude r. The formulation due 
to Pines is free of singularities (except at the origin) as the direction 
cosines are always well defined for any direction in space. Any formulation for 
the potential function and its derivatives is expressed in terms of an infinite 
series . However , it has been found that the gravity model obtained using the 
Pines formulation requires one-third the storage required for the classical 
formulation and that it executes at least 11 times faster due to the eliminatior 
of all trigonometric functions. 

Earth gravity is modeled by the onboard software for all phases of a Space 
Shuttle flight. The degree and order of the gravity model used are varied 
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according to the flight phase and the availability of IMU velocity data. During 
the ascent, entry, and RTLS phases, a fourth degree by fourth order C^x4) 

gravity model is used during coasting periods (IMU delta velocities less than a 

given threshold); otherwise, for these phases, a 2x0 gravity model is used. 
During the onorbit and onorbit checkout phases a 2x0 gravity model is used 
during powered periods (IMU delta velocities greater than a given threshold); 
otherwise, for these phases, a 4x4 gravity model is used. 

B.2.2 ENVIRONMENTAL MODEL 

a. Assumptions - in propagating ' the trajectory, the Shuttle vehicle is treated 

as if it were a point mass as far as gravitational forces are concerned . 

b. Equations - the expansion used to model the Earth's potential V is given by 


V = Z Pn ^ 
n=0 m=0 


where 


n is the degree 
m is the order 

Pn = P ae"r = ag Pn-i/r 

y = Earth's gravitation constant 

Sg = Earth's equatorial radius 

Arj,n,(u) = (n+m-1) (,u) + uAn_i^m(u) 

^1,0 = OJ = 1 

^n,m “ *^n,m^m^®>^^ ^n,m^m^®>^^ 

*^n,m ^n,m classical spherical harmonic coefficients 
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R|jj(s,t) and are the. real and imaginary parts, respectively, 

of the complex number (s + it)“. 

where i = *^'-1 , and s = x/r, t = y/r, u = z/r. 

x,y,z are the Earth- fixed Cartesian coordinates of the Shuttle vehicle 
r is the distance from the center of the Earth to the Shuttle vehicle. 
Rjjj(s,t) and Ijjj(s,t) are given by the recursive relations' 

Rjjj(s,t) = s Rjjj_t (s,t) - t Ij„_-|(s,t) 

I^(s,t) = s IjQ_T(s,t) + t Rjj,_^(s,t) 

where RQ(s,t) = 1, and lQ(s,t) = 0. ' 

The gravitation acceleration of the Shuttle about the Earth is then given by the 
gradient of the potential function V. It can be shown that is given by 

m=0 

- si + t j + uk)) 

where 


CO Q 

Vv = Z 

m=0 a 


n+1 

e 


i,j,k is a set of orthogonal unit base vectors along the coordinate 
axes x,y,z of the Earth- fixed system 

®n,m^®’^^ = ^n,m^m-1^®»^^ ^n,m^m-1 

^n,m^®’^^ “ ^n,m^m-1^®»^^ “ ‘^n,m^m-1 * 

Table B.2-I presents the recommended input gravity model coefficients that are 
baselined for the Space Shuttle software (refs. B.2-3, B.2-4, and B.2-5). 
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TABLE B.2-I GRAVITY MODEL COEFFICIENTS 


Zonal harmonic coefficients 


C(1,0) 

= 

0.0 

C(5,0) 

= 

+2.1923000-07 

C(2,0) 

- 

-1.0826283-03 

C(6,0) 

= 

-5.2317000-07 

0(3,0) 

= 

+2.5418000-06 

C(7,0) 

= 

+3.7219000-07 

0(4,0) 

= 

+1.6086000-06 

C(8,0) 

= 

+1.7564000-07 

Tesseral 

and 

sectorial coefficients 




0(2,1) 


-1.1619-09 

S(2,1) 

= 

-4.1312-09 

0(2,2) 

= 

1 .5654-06 

S(2,2) 

= 

-8.9614-07 

0(3,1) 

= 

2.1625-06 

3(3,1) 

= 

2.6809-07 

0(3,2) 

= 

3.1875-07 

3(3,2) 


-2.1556-07 

0(3,3) 

= 

9-7078-08 

3(3,3) 

- 

1.9885-07 

0(4,1) 

= 

-5.1257-07 

3(4,1) 

= 

-4.4095-07 

0(4,2) 

= 

7.7390-08 

— 3(4,2) 

- 

1.4970-07 

0(4,3) 


5.7700-08 

3(4,3) 

= 

-1.2389-08 

0(4,4) 

= 

-3.4567-09 

3(4,4) 

= 

6.4464-09 

0(5,1) 

= 

-5.8574-08 

3(5,1) 

= 

-7.2105-08 

0(5,2) 

= 

1.0764-07 

3(5,2) 

.= 

5.0363-08 

0(5,3) 

= 

-1.5381-08 

3(5,3) 

- 

-6.4318-09 

0(5,4) 

= 

-1 .9349-09 

3(5,4) 

= 

2.8031-10 

0(5,5) 

= 

4.5428-10 

3(5,5) 

= 

-1 .7529-09 

0(6,1) 

= 

-5.7751-08 

3(6,1) 

= 

-1.1802-08 

0(6,2) 

- 

7.991-09 

3(6,2) 

= 

-4.6527-08 

0(6,3) 

— 

2.3844-10 

3(6,3) 


2.0319-10 
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TABLE Continued 


Tesseral and sectorial coefficients 


C(6,4) 

= 

-3.2820-12 

S(6,4) 

- 

-1.7621-09 

C(6,5) 

= 

-2.2170-10 

S(6,5) 

= 

-4.4098-10 

C(6,6) 

= 

4.0305-12 

S(6,6) 

= 

-6.1204-11 

C(7,1) 

= 

1 .8305-07 

SC7,1) 

= 

1.0137-07 

C(7,2) 

= 

3-4492-08 

S(7,2) 

= 

8.7152-08 

C(7,3) 

= 

2.8003-09 

S(7,3) 

= 

-2.5974-09 

C(7,4) 

- 

-5.9607-10 

S(7,4) 

= 

-2.9899-10 

C(7,5) 


9.3789-12 

S(7,5) 

= 

8.0694-12 

C(7,6) 

= 

-2.1337-11 

3(7,6) 

= 

8.4194-12 

C(7,7) 

= 

1.1576-12 

3(7,7) 

= 

8.9043-14 

C(8,1) 


7.0093-09 

3(8,1) 


3.9788-08 

C(8,2) 

= 

5.0102-09 

3(8,2) 

= 

7.0635-09 

C(8,3) 

= 

-3.8216-10 

3(8,3) 

= 

-6.7434-10 

C(8,4) 

= 

-3.0163-10 

3(8,4) 

= 

3. 7068- 11 

C(8,5) 

= 

-1.0317-11 

3(8,5) 

= 

1.1258-11 

C(8,6) 

= 

-2.6448-12 

3(8,6) 

= 

7.0604-12 

C(8,7) 

= 

3.3552-13 

3(8,7)' 

= 

4.1761-13 

C(8,8) 

= 

-1 .0606-13 

3(8,8) 

= 

8.9361-14 


Earth gravitational radius 
Sq = 2.0925705+07 ft. 

Earth gravitational constant 

U= 1.407645794+16 ft^/sec^ 
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B.2.3 ERROR SOURCES 

As mentioned, the proposed formulation for modeling the gravity field of the 
Earth involves the use of an infinite series. Hence, one source of error is due 
to the truncation of the rather slowly convergent infinite series, the trade-off 
being numerical accuracy versus computation time and storage. Another source of 
error involves the numeric values of the coefficients used in the geopotential 
expansion. The values are based on a batch of observation data and, hence, must 
be consistent with the Earth gravitational constant and radius used in their 
calculation. Many of the higher-ordered coefficients are known to only a few 
significant digits. 
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APPENDIX B.3 
ATMOSPHERE MODELS 


B . 3 . 1 GENERAL DESCRIPTION 

The thermodynamic properties of the real-world atmosphere (temperature, 
pressure, density, and speed of sound) vary widely as functions of altitude, 
latitude, local weather conditions, time of day and year, and solar and 
geomagnetic activity. The following paragraphs describe the various types of 
atmospheric density models used for Space Shuttle development and operational 
support . 


B.3. 1.1 Standard Atmosphere Models 

The properties of the Earth’s atmosphere are described in references B.3-1, 
B.3-2, and B.3-3. Reference B.3-1 is the U.S. Standard Atmosphere, 1962, which 
is the result of a joint effort by NASA, the U.S. Air Force, and the U.S. 

Weather Bureau to depict idealized, middle latitude, year-round mean atmospheric 
conditions. Reference B.3-2, the U.S. -Standard Atmosphere Supplements, 1966, 
was prepared in response to a need for atmospheric tables depicting conditions 
other than the midlatitude mean represented by the tables of the U.S. Standard 
Atmosphere, 1962. In September 1971 the Committee on Extension to the Standard 
Atmosphere (COESA) Working Group reviewed the temperature and density data 
derived from recent satellite and rocket observations. This review revealed a 
need to revise the 1962 Standard at altitudes above 50 kilometers. Data 
available for levels below 50 km were found to be in reasonably good agreement 
with the 1962 Standard. The U.S. Standard Atmosphere, 1976 (ref. B.3-3), with 
tables and graphs extending to 1000 km, was adopted by COESA in February 1975. 
These documents define the various standard atmosphere models, including 
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variation with latitude and season . Also given are- tables of the atmospheric 
properties as a function of altitude, the equations and physical data from 
which the tables were generated, and polynomial approximations to the tabulated 
values . 

B.3.1-2 Jacchia Atmosphere Model 

Reference B.3-4 describes the Modified 1970 Jacchia Atmosphere Model, which 
replaces the 1968 and earlier Jacchia atmosphere models. The primary weakness 
of the 1968 and earlier models was the assumed constant boundary conditions at 
an altitude of 120 km. The 1970 Jacchia model assumes constant boundary 
conditions at 90 km, which closely corresponds to an observed layer of minimum 
variation in the global density distribution.' The following characteristics of 
density are included in this model (ref. B.3-5): 

a. Variation with solar activity (daily and 3 month average) 

b. Diurnal variation 

c. Variation with geomagnetic activity (yearly mean and 3 hourly index) 

d. Semiannual variation 

e. Variation with height 

f. Seasonal-latitudinal variations of the lower thermosphere 

g. Seasonal-latitudinal variation of helium 

To compute density at any point in the atmosphere, the Mission Control Center 
orbit prediction program will employ a model developed by Jacchia (ref. B.3-5) 
that corresponds to the 1970 Jacchia model in analytical structure but whose 
analytical formulation for the geomagnetic effect, semiannual, and helium 
variations have been updated. In the Jacchia model the first four variations 
are computed as changes in the exospheric temperature Too. One may then obtain 
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a rather complicated solution to the diffuse equilibrium equation, knowing T<» 
and the height h. The density p may be determined from the solution of this 
equation for the different gas constituents. Also, the Jacchia model computes 
the last two variations, vrfiich becomes important for altitudes below 300 km. 
Although the model is accurate if requires a large amount of computer storage 
and execution time, thus making it ill-suited for the Shuttle’s onboard computer 
program. However, the onboard program’s atmospheric density model must be 
compatible with the Jacchia model that is used in the ground orbit prediction 
program. 


B.3.1.3 B-M Density Model 

The model described in section B.3«1*2 will be used by the Mission Control 
Center orbit prediction program. This model is quite accurate but it requires 
a large amount of storage and execution time. These undesirable features make 
it ill-suited for the Shuttle onboard program. It is desirable to have an 
onboard model with low storage and time requirements, but mission constraints 
require close agreement between onboard and ground programs. The Babb-Mueller 
(B-M) model was developed with this in mind. 

The model employs a procedure by which coefficients in the model are calibrated 
so that it simulates the Jacchia model. The calibration is accomplished 
premission by direct use of the ground Jacchia model in the I-load program. 

The B-M model is sufficiently general that it can be used with no modifications 
to simulate any reference density model. Additional information about the 
model may be found in reference B.3-6. 
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B.3.2 ONBOARD FORMULATIONS FOR THE SPACE SHUTTLE 

It is anticipated that calibrated models of atmospheric properties will be 
required in three, cases 

a. Surface to 32 km (approx.) - An algorithm to convert sensed atmospheric 
pressure to an equivalent measurement of altitude is required to support the 
barometric altimeter (sec. 4.0). 

b. 30 km to 80 km (approx.) - An algorithm to compute density as a function of 
latitude and altitude will be required if the model presented in reference 
B.3-7 is used for entry analysis as discussed in appendix B.4. 

c. 125 km to 500 km (approx.) - An algorithm to compute density at orbital 
altitudes will be required to support state vector estimation. 

B.3.2.1 Surface to 32 km 

Measurements of the ambient atmospheric pressure by the barometric altimeter are 
converted to measurements of altitude for navigation updating. Section 4.0 
discusses the accuracy with which this may be done . An algorithm for this 
problem is proposed in section 3.5.2 of reference B.3~8. In this algorithm the 
atmosphere is divided into several layers. In each layer the altitude is 
related to the pressure by a polynomial in the logarithm of the pressure ratio 
(P/Pb)- 


N 

h = hb - Hp^b 2 an (ln(p/pb)}^ 


where 

Pb 


base pressure of the layer 
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a^j = polynomial coefficient 
Hp^jj = base pressure scale height 

hp = geopotential altitude from the base of the layer 
p = measured pressure 

The number of layers, and the parameters hjj, Hp, b, a^^, and pj^ .for each lay 
will be provided to the Shuttle. These are based upon ground reduction of 
meterological observations . 

For modeling barometric altimeter measurements in navigation simulations, it is 
both convenient and realistic to ignore the pressure algorithm and to model 
error sources directly as functions of altitude. (See sec. 4.0 for the error 
model and representative numerical values.) 

B . 3 . 2 . 2 30 kilometers to 80 kilometers 

The altitude updating scheme discussed in reference B.3-7 requires an algorithm 
for computation of density as a function of latitude and altitude in the drag 
altitude region 30 km to 80 km. In this model the density is given by a four- 
layer, exponential function of altitude. 

P = Poe-^/\ 

where 

Pq = base density 
hg = scale height 

h = altitude above the Fischer’s ellipsoid 

In the above model, the base density (p^) and scale height (hg) are constants 
for each layer. Reference B.3-9 gives the following base density and scale 


B-17 



ORlGSi^AL PAUE 6^ 

OF POOR QUALITY 

79FM5 

height constants used in modeling a four-layer density model for navigation 
error studies. 


FOUR LAYER DRAG FORMULATION 


Altitude region, feet 

Base density, kg/m^ 

Scale heights, meters 

85 200 to 155 348 

1.6546565 

6705.8 

155 348 to 172 010 

.47779728 

8135.6 

172 010 to 202 069 

.42736363 

8278.9 

202 069 to higher 

1.1927824 

7274.9 


However,- it is possible to load a seasonal atmosphere where the base density 
and scale height are constants for each layer, but vary in value depending on th 
season and latitude of interest (ref. B.3-7). The following atmospheres have 
been specified for use 

Latitude Season 

1962 Standard Atmosphere Model All 

60° N Winter (January) 

Summer (July) 

30° N Winter (January) 

Summer (July) 

Spring (April) 

These models were obtained from reference B.3-1 and exhibit an average curve fit 
error of 3 percent, with the maximum curve fits error never greater than 8 
percent. Tables B.3-1 through B.3-IV give the altitude regions, densities, base 
densities, and scale heights for each density model referenced above. The 
errors associated with the density model proposed for this altitude region 
can be grouped into two categories . 
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a. Diurnal variation, which is a variation during a single day due to the 
heating of the atmosphere 

b. Variations from the model during any month of that season 

These errors can be incorporated by varying Pq and hg or by specifying Ap/p 
to account for the known excursions from the density model. The following 
tables present some recommended three-sigma error values for these error sources 

Diurnal Variation ^0 kilometers 

All latitudes Ap^ = 2.756;Ahg = .05? 

Variation from seasonal (for altitudes from 30 kilometers to 80 kilometers. 


Latitude 

Season 

APp/Pp, percent 

Ah„/h„, percent 

30° N 

Winter 

(January) 

12.01 

1.39 


Summer 

(July) 

9.90 

.87 


Spring 

(April) 

34.43 

10.91 

60° N 

Winter 

(January) 

4.32 

2.31 


Summer 

(July) 

7.92 

1.68 

1962 Standard 

Atmosphere Model - 

All . 10.7 

4.69 


Table B.3-V presents the three-sigma density percentages versus altitude (30 to 
80km) for the reference atmospheres. These two error sources are Independent so 
that during any season the diurnal variations exist along with the predicted 
variation for that season. Also, it can be seen that the use of a seasonal 
atmosphere instead of the 1962 Standard Atmosphere reduces the predicted error 
(ref. B.3-7). 


65 kilometers 

APq =5.34?;Ahg = .10$ 
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B-3.2.3 125 kilometers to 500 kilometers 

In order to predict trajectories of orbiting vehicles, a model of the 
atmospheric density is required in the computation of aerodynamic forces. 
Reference B.3-10 discusses the proposed replacement of the baselined Russian 
Atmospheric density model with the B~M density model for simulating the ground 
Jacchia model in the onboard software. The B-M density model, which computes 
density at orbital altitudes, uses the following equation for modeling density 
(ref. B.3-5). 

p = p^exp (A+B) 

where A is the nighttime vertical profile 
A = a^^ + a 2 h + ^ 3 /^ 
and B is the diurnal effect 

B = (b-j+b 2 h+b 3 /h) ((1 + cos 

where 

pQ = base density 
h = altitude 

a^ ,a 2 > a^jb-j ,b 2 ,b 3 = model parameters 

M = power exponent of the diurnal term 

cos\|; = (z sindg + cosfigCXcosy + Ysiny)]/R 

R = vehicle's distance from the center of the Earth 

6 g = declination of the Sun 

y = right ascension of the diurnal bulge 

X,Y,Z = coordinates of the vehicle in M50 coordinates 
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The seasonal- latitudinal variation may be included in the B-M model. Since the 
Jacchia model computes this term as an explicit variation in the density, one ma 
incorporate this in the B-M density model. The seasonal-latitudinal variation 
that is given by Jacchia is (ref. B.3-5) 

p = p*.10e 

where 

e = 0,02 (h-90) (6^/1 6yl ) exp (-0.45(h-90)] sin^^v x 
sin (360 (d+100)/y) 

where 

h = altitude in kilometers 
d = number of days into a year 

Y = number of days in a year 

- declination of the vehicle 

This term may be rewritten so that the final form of the model is 
p = PqSxp (A+B+C) 

where 

C = 0.046051 Z| (Z/R2)(h-90)exp(-0.i}5(h-90)] sin(360(d+100)/Y) 

where 

Z = Z component of the vehicle's position vector 
h = altitude in kilometers 

Y z days in the year 
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d = initial day in the year 

R = vehicle’s distance from the center of the Earth 

The sine term may be computed once with the initial day (d) and days in the 
year (Y) and then assumed constant. Reference B.3-5 discusses the calibration 
procedure for determining the nighttime vertical profile and diurnal coeffic- 
ients through calibration to the Jacchia model. If one assumes that the power 
exponent of the diurnal term (M) , Pq and the bulge angle <i) are known 
constants and unaffected by the solar activity or geomagnetic changes, the B-M 
density model has six coefficients to be determined. By careful study of the 
Jacchia model a value of <{) = 37** > M = 2.75, and Pq = 1.22499 kg/m has been 
adopted (ref. B.3-5). 

However, all these values may be refined to obtain closer agreement with the 

Jacchia model. To compare the B-M density model to the Jacchia model, reference 

B.3-5 conducted several numerical experiments. Instead of comparing directly 

computed densities , predicted satellite positions are compared from a numerical 

orbit computation program. Several different orbits were chosen for the 

comparison; thus, the position difference is a good indication of the global 

difference between the density models used in the experiments. For each orbit- 

2 

the ballistic number was an average value for the Shuttle, BN = 100 Ib/ft , 
and the coefficient of drag, was set to 2.2. Table B.3-VI displays the 

chosen orbits and table B.3-VII presents the results. The results in table 
B.3-VII are the differences in predicted position using the Jacchia model as 
compared to using the other three models or neglecting drag completely. It is 
apparent from the results of table B.3-VII that the B-M model gives better 
agreement to the Jacchia than either the U.S.S.R. or ADMB models. 
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TABLE B.3-I. DENSITY MODEL - 


Altitude, 

m 

P’3 

kg/m-^ 


30480.0 

0.01710 


32308.8 

.01292 


34137-6 

.00940 


35966.4 

7.295 X 

10-3 

37795.2 

5.533 X 

10-3 

39624.0 

4.221 X 

10-3 

41452.8 

3.239 X 

10-3 

43281.6 

2.498 X 

10-3 

45110.4 

1.937 X 

10-3 

46939.2 

1.439 X 

10-3 

48768.0 

1.142 X 

10-3 

50596.8 

9.066 X 

10-4 

52435.6 

' 7.219 X 

10-4 

54254.4 

5.796 X 

10-4 

56083.2 

4.640 X 

10-4 

57912.0 

3.702 X 

10-4 

59740.8 

2.945 X 

10-4 

61569.6 

2.345 X 

10-4 

63398.4 

1.875 X 

10-4 

65227.2 

1.489 X 

10-4 

67056.0 

1.174 X 

10-4 

68884.8 

9.186 X 

10-5 

70713.6 

7.128 X 

10-5 

72542.4 

5.482 X 

10-5 

73761.6 

4.578 X 

10-5 

75590.4 

3.464 X 

10-5 

77419.2 

2.593 X 

10-5 

79248.0 

1.909 X 

10-5 

81076.8 

1.351 X 

10-5 

82905.6 

9.561 X 

10-6 

84734.4 

6.766 X 

10-6 


OF 




B-23 


79FM5 



ORIGINAL PAGE IS 
OF POOR QUALITY' 


79FM5 


TABLE SEASONAL DENSITY MODEL - JANUARY 60® N AND JULY 60° N 


JANUARY 60® N 


JULY 60® N 


Alti- 

tude, 

km 

kg/m- 

1 

kg/m3 

hg, 

m 

Alti- 

tude, 

km 

P> . 

kg/m- 

5 

kg/ra^ 

hg, 

m 

32 

1.170 

X 

10-2 



32 

1.443 

X 

10-2 



34 

8.487 

X 

10-3 



34 

1.061 

X 

10-2 



36 

6.106 

X 

10-3 

1.8898 

6301.1 

36 

7.857 

X 

10-3 

1.525 

6881.0 

38 

4.424 

X 

10-3 



38 

■5.861 

X 

10-3 



40 

3.228 

X 

10-3 



40 

4.402 

X 

10-3 



42 

2.371 

X 

10-3 



42 

3.328 

X 

10-3 



44 

1.753 

X 

10-3 



44 

2.548 

X 

10-3 



46 

1.303 

X 

10-3 



46 

1.969 

X 

10-3 



48 

9.571 

X 

10-4 



48 

1.525 

X 

10-3 



50 

7.336 

X 

10-4 

.7351 

7286.3 

50 

1.192 

X 

10-3 

.5824 

8100.4 

52 

5.641 

X 

10-4 



52 

9.315 

X 

10-4 



54 

4.338 

X 

10-4 



54 

7.330 

X 

10-4 



56 

3.377 

X 

10-4 



56 

5.792 

X 

10-4 



58 

2.619 

X 

10-4 



58 

4.561 

X 

10-4 



60 

2.014 

X 

10-4 



60 

3.613 

X 

10-4 



62 

1.539 

X 

10-4 

.5711 

7533.7 

62 

2.869 

X 

10-4 

.4554 

8399.2 

64 

1.175 

X 

10-4 



64 

2.259 

X 

10-4 



66 

8.955 

X 

10-5 



66 

1.761 

X 

10-4 



68 

6.820 

X 

10-5 



68 

1.359 

X 

10-4 



70 

5.221 

X 

10-5 



70 

1.038 

X 

10-4 



72 

4.008 

X 

10-5 



72 — 

7.824 

X 

10-5 



74 

3.062 

X 

10-5 

.7367 

7321 .0 

74 

5.821 

X 

10-5 

3.3202 

6704.3 

76 

2.327 

X 

10-5 



76 

4.268 

X 

10-5 



78 

1.759 

X 

10-5 



78 

3.080 

X 

10-5 



80 

1.323 

X 

10-5 



70 

2.182 

X 

10-5 
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TABLE B.3-III.- SEASONAL DENSITY MODEL - JANUARY 30® N AND JULY 30 
JANUARY 60® N JULY 60® N 


Alti- 
tude , 
km 

kg/m^ 

kg/ m3 

m 

Alti- 

tude, 

km 

P, 

kg/m3 

30 

1.730 

X 

10-2 



30 

1.856 

X 

10-2 

32 

1.265 

X 

10-2 



32 

1 .362 

X 

10-2 

34 

9.273 

X 

10-3 



34 

1 .003 

X 

10-2 

36 

6.840 

X 

10-3 

1.5164 

6706.4 

36 

7.425 

X 

10-3 

38 

5.076 

X 

10-3 



38 

5.530 

X 

10-3- 

40 

3.788 

X 

10-3 



40 

4.141 

X 

10-3 

42 

2.843 

X 

10-3 



42 

3.118 

X 

10-3 

44 

2.145 

X 

10-3 



44 

2.360 

X 

10-3 

46 

1.627 

X 

10-3 



46 

1.795 

X 

10-3 

48 

1.250 

X 

10-3 



48 

1.384 

X 

10-3 

50 

9.697 

X 

10-4 

.5853 

7844.6 

50 

1.070 

X 

10-3 

52 

7.576 

X 

10-4 



52 

8.433 

X 

10-4 

54 

5.944 

X 

10-4 



54 

6.634 

X 

10-4 

56 

4.646 

X 

10-4 



56 

5.20 

X 

10-4 

58 

3.618 

X 

10-4 



58 

4.061 

X 

10-4 

60 

2.818 

X 

10-4 



60 

3.180 

X 

10-4 

62 

2.191 

X 

10-4 

.5841 

7846.8 

62 

2.491 

X 

10-4 

64 

1.693 

X 

10-4 



64 

1.936 

X 

10-4 

66 . 

1.299 

X 

10-4 



66 

1.492 

X 

10-4 

68 

9.893 

X 

10-5 



68 

1.140 

X 

10-4 

70 

7.479 

X 

10-5 



70 

8.632 

X 

10-5 

72 

5.609 

X 

10-5 



72 

6.469 

X 

10-5 

74 

4.171 

X 

10-5 

2.5412 

6679.2 

74 

4.796 

X 

10-5 

76 

3.074 

X 

10-5 



76 

3. .51 4 

X 

10-5 

78 

2.244 

X 

10-5 



78 

2.543 

X 

10-5> 

80 

1.597 

X 

10-5 



80 

1.815 

X 

10-5 



1.5412 


.6048 


.5656 


3.0678 


79FM5 
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6778.4 


7933.4 

8009.3 


6645.7 
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TABLE SEASONAL DENSITY MODEL - APRIL 30® N 


Altitude, ' p, 
km kg/m^ 


^ “s' 

kg/m^ m 


30 

1.74 

X 

10-1 



35 

8.22 

X 

10-3 

1.2536 

7013.8 

ho 

4.02 

X 

10-3 



45 

2.05 

X 

10-3 



50 

1.10 

X 

10-3 

.4794 

8249.9 

55 

6.10 

X 

10-4 



60 

3.40 

X 

10-4 

.5176 

8156.1 

65 

1.79 

X 

10-4 



70 . 

9.03 

X 

10-5 



75 

4.35 

X 

10-5 

2.6028 

6781.6 

80 

1.96 

X 

10-5 




OF POOR QOAUT'? 
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TABLE B.3-V.- 3cr A6 




VERSUS ALTITUDE FOR ENCLOSED ATMOSPHERES 


Altitude 

(km) 

1962 

Standard 

atmosphere, 

percent 

60° N 
January, 
percent 

60° N 
July, 
percent 

30° N 
January, 
percent 

30° N 

July, 

percent 

30° N 
April , 
•percent 

30 

10.09 

4.83 

2.13 

7.03 

7.38 

5.65 

32 

11.61 

5.42 

1.71 

6.69 

7.15 

3.09 

34 

13.10 

6.01 

1.30 

6.34 

6.91 

.60 

36 

14.57 

6.60 

.88 

6.00 

6.67 

1.83 

38 

16.01 

7.18 

.47 

5.65 

6.44 

4.21 

40 

17.42 

7.76 

.06 

5.31 

6.20 

6.52 

42 

18.82 

8.33 

.35 

4.96 

5.97 

8,78 

44 

20.18 

8.90 

.76 

4.61 

5.74 

10.99 

46 

21.53 

9.47 

1.16 

4.26 

5.50 

13.14 

48 

22.85 

10.03 

1.57 

3.90 

5.27 

15.24 

50 

24.16 

10.59 

1.97 

3.55 

5.04 

17.29 

52 

25.43 

11.15 

2.37 

3.20 

4.81 

19.29 

54 

26.69 

11.70 

2.77 

2.84 

4.58 

21.24 

56 

27.93 

12.25 

3.17 

2.48 

4.35 

23.15 

58 

29.14 

12.80 

3.56 

2.12 

4.12 

25.00 

60 

30.34 

13.34 

3.96 

1.76 

3.89 

26.82 

62 

31.51 

13.88 

. 4.35 

1.40 

3.66 

28.59 

64 

32.67 

14.42 

4.74 

1.04 

3.43 

30.31 

66 

33.81 

14.95 

5.13 

.67 

3.20 

32.00 

68 

34.92 

15.48 

5.52 

.31 

2.97 

33.64 

70 

36.02 

16.01 

5.90 

.06 

2.75 

35.25 

72 

37.10 

16.53 

6.29 

.43 

2.52 

36.81 

74 

38.16 

17.05 

6.67 

.80 

2.30 

38.34 

76 

39.20 

17.57 

7.05 

1.17 

2.07 

39.83 

78 

40.23 

18.08 

7.43 

1.54 

1.85 

41.29 

80 

41 .24 

18.59 

7.81 

1.91 

1.62 

42.71 
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■ TABLE B.-3- 

-VI.- ORBITS 

USED IN 

PREDICTION 

EXPERIMENTS 

■ 

- 

A 

■ B 

C 

- D 

E 

perigee (km) 

220 

300 

300 

220 

Same as 

apogee (km) 

380 

600 

600 

380 

A except 

eccentricity 

012 

.022 

.022 

.012 

epoch 

period (min) 

90.5 

93.6 

93.6 

90.5 

is 

argument of perigee 

0 

0 

180® 

0 

12:00 

ascending node 

0 

0 

0 

0 

January 1 , 

inclination 

30 ® 

30 ® 

30 ® 

90® 

1977 

epoch 

12:00 


January 

1, 1975 
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TABLE B.3-VII.- POSITION DEPENDENCE ON DENSITY MODEL 


Orbit time of integration, Position difference, km 

days 




No drag 

U.S.S.R. 

AMDB 

B-M 

A 

0.5 

8.9 

0.2 

0.3 

0.005 


1.0 

36.1 

.8 

1.5 

.018 

B 

0.5 

.83 

.03 

.07 

.004 


1.0 

3.16 

- .07 

.28 

.015 

C 

0.5 

.5 

<.001 

.01 

.004 


1.0 

1.8 

.02 

.04 

.018 

D 

0.5 

8.6 

.1 

.3 

.008 


1.0 

34.4 

.5 

1.5 

.027 

E 

0.5 

13.9 

3.4 

0.4 

.071 


1.0 

56.0 

19.2 

1.9 

.275 


B-29 




79FM5 


B . 3 . 3 Re ference s 

B-3-1 U.S. Committee on Extension to the Standard Atmosphere: U.S. Standard 

Atmosphere, 1962. Government Printing Office (Wash., D.C.), 1962. 

B.3-2 U.S. Committee on Extension to the Standard Atmosphere: U.S. Standard 

Atmosphere Supplements, 1966. Government Printing Office, (Wash., D.C.) 
1962. 

B.3-3 U.S. Committee on Extension to the Standard Atmosphere: U.S. Standard 

Atmosphere, 1976. Government Printing Office (Wash., D.C.), 1962. 

B.3-4 Lewis, J. R-; Sabayrac, R. B. ; Austin, G. A.; Cockrell, B. F.: RTCC and 

MOPS Requirements for Modified 1970 Jacchia Atmosphere Model. MSC 
IN 70-FM-202, January 26, 1971. 

B.3-5 Mueller, A. C. : Atmosphere Density Models. Analytical and Computa- 

tional Mathematics, Inc., June 1977. 

B.3-6 Mueller, A. C.: Atmospheric Density Models. Analytical and Computa- 

tional Mathematics, Inc., January 25, 1977- 

B.3-7 Purifoy, D. D.: New Wind Bias, Density anbd Density Error Models for 

OFT-1 Simulation. JSC Memorandum FM83 (77-337), September 6, 1977- 

B.3-8 Clifford, J. B., Jr.; and Cockrell, B. F.: Space Shuttle Navigation 

Software , Volume II - Environment Routines . MSC IN 72-FM-132 , June 1 , 
1972. • 


B-30 



79FM5 


B.3-9 McClain, C. R.: Comparison of Two-IMU Nominal Navigation Performance - 

Between the Single and Four-Layer Atmospheric Density Models, MDTSCO 
Working Paper No. E914-8A/B-014, November 5, 1976. 

B.3-10 Lukaszewski, D. A.: Onorbit Navigation Software Atmospheric Density 

Model {Proposed Replacement for Baselined Russian Design), MDTSCO 
Transmittal Memo No, 1.4-MPB-80l|, August 30, 1977. 


B-31 



79FM5 


APPENDIX B.4 

DRAG ALTITUDE PSEUODOMEASUREMENT 
B.4.1 GENERAL DESCRIPTION 

Vertical channel instability is an. inherent problem with all unaided inertial 
navigation systems. The Shuttle navigation system experiences this problem 
during atmospheric reentry. External navigation aids, such as TACAN, are 
generally not available during the initial phase of entry due to blackout 
ionization of the air. 

During entry only IMU data are available. Vehicle aerodynamic characteristics, 
an atmospheric density model, and IMU contact acceleration are used to compute 
altitude for processing with a Kalman filter to stabilize the vertical 
navigation channel. The measurement is referred to as pseudo since it is not 
directly measured. All the acceleration measured by the IMU accelerometers 
during entry is assumed to be attributed to aerodynamic forces. The following 
sections describe the drag altitude pseudomeasurement simulation approach and 
error model. 
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B.H.2 SIMULATION APPROACH 

On nominal trajectories appreciable atmospheric drag is encountered at 
altitudes below 260 000 feet. Drag measurements begin vrtien the accelerometers 

P 

sense 11 ft/ sec and continue until the aerodynamic probe is extended and 
barometric altimeter data can be processed. This occurs at an altitude of 

approximately 100 000 feet. The drag updating technique is an effective method 

for reducing the error buildup during the blackout portion of entry. During the 
drag update period the altitude pseudomeasurement is calculated from 
accelerometer outputs, a model of air density as a function of altitude, and a 

model of drag as a function of the angle of attack. Figure B-1 presents a 

block diagram showing the computation of the drag altitude pseudomeasurement. A 
functional description of each block is as follows. 

a. Block 1 - computes the angle of attack (a) to be used in calculating the drag 
coefficient by the following equations 

Xr = Vr > Ep (Ip X V) 

ZrB = TbM50^ 1 b 
a = arctan (Vpp (3)/Vp3(1)) 


where 


V = the Shuttle’s velocity vector in M50 coordinates 
Ep = Earth rate 

Ip = unit vector along the Earth's pole 

“ transformation matrix from body to M50 coordinates 
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Typical values for Earth rate and Earth's pole unit vector (ref. B.4-1) are 

= 0.72921151 E-4 rad/sec 

0.28324166 E-2 
Ip = -0.50463837 E-4 
0.99999598 

b. Block 2 - computes the drag coefficient (C^^) as a function of the angle of 
attack and premission determined coefficients (Cq, C^, C^) by the 
following equation: 

Cjj = Cq + C^a + 

where 


a = flight angle of attack in radians 

Cq, C^, = experimentally determined coefficients from wind 

tunnel tests for calculating the drag coefficients 

The following typical values for calculating the drag coefficient are given in 
reference B.4-1. 

Cq = 0.2452220322 

= -0.831161 rad"^ 

C2 = 2.4092462 rad"^ 

c. Block 3 “ calculates the ballistic coefficient (3) 

3 = CpA/W 
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where 

A = frontal area of the Shuttle 
W = mass 

= drag coefficient (unitless) 

Due to the longitudinal center of gravity constraint the Orbiter's weight at 
entry, will be kept constant at approximately 183 404 pounds (ref. B.4-2). 

d. Block 4 - computes the density (p) by the following equation. 

P = 2|lyp • /(Vr23At) 

where 

Vr/Vr 

the accelerometer sensed change in the velocity vector over, the 
time interval AT 

time interval between accelerometer samples 
velocity vector of Shuttle' relative to the atmosphere 

|v„| 

calculated in block 3 

e. Block 5 *- computes the drag pseudoaltitude as a function of density by the 
following equation. 

h = hg ln(pQ/p) 

where 


-VR 

AV 


AT = 


^R 


8 = 


hg = atmospheric density scale height 
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Pq = base density 
p = density of the atmosphere 

Section B. 3 . 2.2 of appendix B.3 discusses the baselined four-layer atmospheric 
density model. In this model the scale heights and base densities vary 
according to estimated altitude. 

B. 4.3 ERROR MODEL 

In the navigation filter the error due to drag measurements is modeled as an 
ECRV. The atmospheric density and drag coefficient model errors comprise the 
bias error, vrtiile the accelerometer and atmospheric turbulence errors make up 
the random error seen in the estimate of drag altitude. The following error 
model (1 -sigma values) presented in reference B.4-3 is used to corrupt the 
altitude calculated in block 5 . 

Drag 
Drag 
Drag 


altitude measurement bias 
altitude measurement random noise 
altitude measurement bias time constant 


3750 

ft 

250 

ft 

720 

sec 
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B.5 REFRACTION CORRECTION MODELS 


B.5.1 GENERAL DESCRIPTION 

A complete refraction model would account for tropospheric and (for orbital 
altitudes) ionospheric refraction effects on the range, range rate, integrated 
Doppler, and elevation angle data for all ranges and elevation angles. These 
effects can be well modeled by a ray tracing program (ref. B.5-1) assuming that 
the physical parameters of the atmosphere are- known . This is done by treating 
the atmosphere as a sequence of layers with the refractivity held constant in 
each layer. The refractivity of a given layer is determined from the single 
exponential refractivity profile. The spherical form of Snell's law is then 
used to trace the ray path through the atmosphere. The chief disadvantage of 
this method is that the computations require a considerable amount of computer 
time and are, therefore, unsuitable for general simulation programs or for 
onboard computer use. In these applications less complex and faster 
computational techniques are required. 

A practical refraction model algorithm is discussed below, with estimates of the 
magnitudes of the effects and errors remaining after corrections have been 
applied. It is generally concluded that no matter how sophisticated the 
refraction correction model is, it should not be used for high-quality 
navigation work at elevation angles below 2 or 3 degrees when the vehicle is at 
large distances from the station. Because of diurnal, hourly, vertical, and 
horizontal variations of the index of refraction it is virtually impossible to 
construct an accurate, longrange , low-angle, refraction correction. 
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B.5.2 TROPOSPHERIC REFRACTION MODELS 

The following discussion of tropospheric refraction is taken directly from 
reference B.5-2. Figure B-2 shows a simplified model of the atmospheric 
refraction of electromagnetic waves. 

Two effects are involved. First, the electromagnetic waves travel more slowly 
through the atmosphere than they do through a vacuum. The ratio of the speeds 
is called the index of refraction, n. Due to the slower velocity in the 
atmosphere it takes longer for waves to travel from the tracking station to the 
vehicle than they would in a vacuum. This makes the measured distance appear 
greater than the actual distance. Secondly, the electromagnetic waves are bent 
at the outer edges of the atmosphere . This makes the path followed by the waves 
longer than the geometric, straight-line distance. This makes the measured 
distance appear to be greater than the actual distance p . 

Before proceeding, let us precisely define a few of the terms shown in figure 
B-2. 


= the index of refraction for the electromagnetic wave at the 

tracking station antenna, ng-1 is called the refraction modulus, 
or refractivity. 

Hg = Atmospheric decay constant, or atmospheric scale height. Looking 
straight up it is the equivalent height of a slab atmosphere whose 
index of refraction is ng in every case. 

H = altitude of the vehicle above the tracking station. 
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Hq = radius of curvature of the Earth in the vicinity of the tracking 
station. 

If Eg = 6 378 166 meters is the equatorial radius, e = 1/298.3 is the ellipticity 
or flattening, of the Earth, and (|) = geodetic latitude; then 


Rc = 


(1 -• e)^ Re 

(l-(1-(1- e)^) sin^)3^^ 


- + ALT 


station 
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From the Rq, Pg^ triangle and the law of sines, we see that 

Rc Rq + Hg 

sin 0Q sin(90° + Ej^) 

or (1) 


1 

sin 6a = cos Em 

^ 1 4. Hg/Rg M 


Also from this same triangle, using the law of cosines, we can show that 


Pa = t Em - sin E^J 


( 2 ) 


or , alternatively 

(2 + Hg/Rg)Hg 

Pa = 

sin Em + (1 + Hg/Rc)^ - cos^Em 

If ng is the index of refraction of the slab atmosphere, then from Snell's 
law of refraction 

sin 0^, = ng sin 0g 

Substituting equation (1) into the above gives 

”s 

1 + Hg/Rg 

Note that for ng > 1 + Hg/Rg it is possible to have reflection rather than 


refraction . 
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From the p^, pjj, p triangle .we see that 

p^ sin (0jj - Gg) = p sinCe^j-e^-(Ej^E) ) (4) 

, 2 

pjj = p<^ + p^ _ 2PPa gos(Em-E) C5) 

This now gives us five nonlinear equations in the five unknowns: 0^, 

6jj, pg, and Pjj. The quantities p and E are the independent variables 
or given in quantities. There is no analytic solution of these equations; 
however, an iterative solution is easily obtained as shown below. 


% = E 


D([) a I = 1,3 


0, = arc sin 


cos Ejj 
1 + Hg/Rc 


= arc sin 


ng cos Ef^ 

1 + Hg/Rg 


Pa = Rc + Hg/Rg)^ - cos^Em - sin Em] 


Pa 

a Ej^ = E + (0jj - 0^) - arc sin -- sin(0j^ - 0^) 


p|3 = p2 + p2 _ 2pp^ cos(Em-E) 

H 


The refraction correction for range. Zip, is measured range, + Pb' 

minus computed range , p . Thus 
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Ap = ngp^ + Pb + P 

The previous refraction correction equations are only good for the vehicle well 

above the Earth’s atmosphere. As the vehicle enters the atmosphere, part of the 

atmosphere will be above the vehicle, leaving a lesser quantity of atmosphere 

below the vehicle. Thus, for a vehicle in the atmosphere we would like to 

lower the scale height, Hg. This can be done by replacing Hg 
* 

with Hs , 

Hg = [ 1 - exp(-H/Hg)]Hs (6) 

where altitude above the tracking station ,is conveniently calculated from 

H = + 2PRq sin E - Rc (7) 

» » 

Note that as H becomes small, Hg H. Except for H = 0, Hg is always 

less than H, for H positive. That is, the vehicle will always be above 

the slab atmosphere for H positive. Also note that as H becomes large, 

* 

Hg Hg. As it turns out, to a first approximation the refraction modulus as a 
function of altitude is given by (ng-1) expC-H/Hg). This being the case, we 
can show that eq. (6) is precisely the right equation to use to modify the scale 
height . 

* 

If we substitute Hg for Hg in our previous equations we encounter two 
difficulties. One is for negative altitudes and the other is the reflection 
problem. We can easily resolve the negative altitude problem by noting that 
when H is negative (station on a hill during landing) the range p will 
never be very large. Thus when H is negative we will use the approximate 
refraction corrections. 
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% ~ E 

for H < 0 

Ap =; (ng-DP 

The other problem concerns eq. (3) > 

ng cos Ej^ 

0K = arc sin - 

^ 1 t HJ/Rc 

* 

As mentioned earlier, for 1 + %/Rq > ng it is possible to have reflection 

« 

rather than refraction. This will occur when Hg ~ 2000 meters and H~2500 
meters. At this altitude we can obtain a reflection from the upper part 
of the slab atmosphere rather than a refraction. While this effect can occur 
in nature, and cause a mirage, it is not something that we desire to model 
here. In fact, it will not occur in an atmosphere whose refraction modulus 
is smoothly varying according to (ng-1)expC-H/Hg) . The reflection is due 
to the abrupt discontinuity in the refraction modulus at the top of the slab 
atmosphere. When reflection occurs, will-be slightly greater than 90 

degrees . We can not have > 90° because in this ease there will be regions 
in the sky where we would not be able to see the Shuttle. To resolve our 
difficulty we will limit ^ by checking 

ng cos 

Whenever this quantity is greater than 1 , we will reset it equal to 1 , causing 6jj 
be computed as 90 degrees. 

Thus, the suggested algorithm for computing the range and elevation angle 
refraction corrections is summarized in figure B-3- 
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H =Vr0 + + 2 pR(;; sin E -Rq 

% = E 

IF (H > 0 ) GO TO a 
Ap = (ng - Dp 
RETURN 


a 


b 


Hs = [r - exp(-H/Hs)]Hs 
DO b I = 1,3 


cos El 


0* = arc sin 

3. 


M 


1 + H^/Rg 


A = 


ng cos 


IF (A > 1) A = 1 


0jj = arcsin (A) 

Pa = -\|(1 + Hs/Rc)^ - cos 2E{^ - sin Ej^] 


Ej^ = E + (0jj - 0^) - arcsin 


fPa 
!in — 

Lp 


sin (0^ - 0_) 
b ^ 


3 


Pb 




Pg - 2pp^ cos (Efj-E) I 


Ap = np^ + pjj - p 


RETURN 


Figure B-3-- Suggested algorithm for computing refraction 
corrections for the Space Shuttle. 
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B.5.3 IONOSPHERIC REFRACTION MODELS 

Ionospheric refraction affects signals between orbiting vehicles and the 
ground. It is a function of the electron density in the ionosphere. The 
electron density is, in turn, a function of the solar and geomagnetic activity 
parameters and the relative positions of the Sun, vehicle, and station. 

A model for the density similar to that for the exospheric temperature (appendix 
B.3) might be available. If mean values for the electron density and the 
ionospheric limits are used instead, the correction for ionospheric refraction 
will not be more accurate than 20 percent. 


The following equation, recommended by Cubic Corporation (refs. B.5-3 and 

may be used to compute the effect of ionospheric refraction on range 
measurements derived from frequency modulation. 


Ap = 


1325 HqMp2 


f^W(.05 cosEq + sinE^) 


[tan“% - tan“^ 


W 






(ft) 


where 


H = the altitude of the vehicle (ft) 

Mp 2 = the maximum electron density in the F 2 ionospheric layer (ft“^) 
f = carrier frequency (Hz) 

W = 

Eq = incident elevation angle 
Hq = height of Mp 2 layer (ft) 

%>% = upper and lower heights of the half values of the electron 
density (ft) ■ • ' 
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Typical values are 

Mp2 = 13.4987 X 10^ 

= 725 000 ft 

Hji = 1 368 600 ft 

Hl = 608 266 ft 

The ionospheric refraction correction is typically about 75 to 80 percent 
of the tropospheric correction for a frequency of about 1.5 GHz. Because 
the correction is inversely proportional to the square of the frequency, 
it becomes small at the higher frequencies in the S-Band. 

b.5.4 refraction error model 

The real atmosphere is, of course, vastly more complicated in its refractivity 
structure than the simple refractivity profile given above. One of the chief 
differences is that the model assumes horizontal homogeneity. That is, every 
place with the same height has the same refractivity. Clearly, this is not 
true in the real world. 

The model for the effects of atmospheric variations on range and Doppler 
refraction corrections is modeled in the Navigation Analysis Program (ref. 
B.5-5) as a scaling error k, which is the sum of a random variable that 

is exponentially correlated over spacecraft position changes,' and is common 

» 

to all stations, and a statistically independent random variable. The modeled 
refraction correction is, therefore, given by the expression 

Ap = (1 + k) Ap 

k = e + q 
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v*.ere 

Ap is the refraction correction computed from the simplified 
atmosphere model 

k is the random scaling error to account for atmospheric variations 
e is an exponentially correlated random variable 

Ti is' white noise 

The bias scale factor £ is computed as a first-order process in Ar. 

^i+1 = exp(-Ar/T)£j^ + a^jj^(l_exp(-2Ar/T))1^2 

inhere 

= |Si4.1 - Si I 

T is the correlation distance constant 

standard deviation for the random variable 

Jumerical values for t and difficult to obtain. A reasonable 

/alue for probably on the order of 5 percent for elevation angles 

ibove 3 degrees and 10 to 20 percent below 3 degrees. 
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B.5-4 
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B.6 LIGHT TIME ALGORITHM* 

The Light Time Algorithm is necessary because the actual times of signal 
reception/transmission at the vehicle ty and transmission from the ground based 
transmitting station t^ are not necessarily known. To calculate the range and 
range rate measurements it is necessary to compute the acutal time required for 
the signal to travel to and from the vehicle. The. times are obtained as 
follows. 

A signal is transmitted from a ground based station at time and travels 
to the vehicle at the speed of light . This signal is received and .transmitted 
by the vehicle at time ty to a ground based receiving station. The onboard 
processing time is assumed to be zero. The input range and range rate obser- 
vation data is tagged only with the time of data reception tj,. In order to 
accurately compute the observations it is necessary to calculate the trans- 
mission times, At 2 = tj, - ty and At-j = ty - t^, and in turn the times ty and t.^. 
The transmission times are obtained from the Light Time Algorithm. This method 
determines the transmission times by successive approximations. It considers 
two transcendental equations 


.pp(t.y) ~ Atp) 

At2 = = — ^ (G1) 


P-|(tt) Pl^ty - At-l) 

At. = = (G2) 

c 


*This appendix is taken directly from the following document: "Houston Operations 

Predictor /Estimator (HOPE) Engineering Manual - Revision 01,” TRW Report 
11176-Hi|08-RO-01, June 1970. 
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where 

P 2 = magnitude of the downlink range vector P 2 
p.j = magnitude of the uplink range vector P-j 
c = speed of light. 

The equation for the downlink range vector P 2 iS' 

-> ->■-»- 
P2 = r - r^ 

where 


r = position vector of the vehicle at time t^ 

= position vector of the ground based receiving station at time 
^r' 

The equation for the uplink range vector p.j is 


Pi = 



(G-4) 


where 


r = position of the vehicle at time t^ 

r^ = position of the ground based transmitting station at time t^. 


Equations (G-1) and (G-2) must be solved by an iterative procedure, since 
they are transcendental equations. Equation (G-1) begins with the known 
time tp and retraces the received signal backwards along the light path to 
solve for t^. Equation (G-2) is then used to retrace the signal from the 
vehicle backwards to the ground based transmitting station to solve for t^. 
The solution to these two equations is obtained as follows. 


4 
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First, the. solution to the downlink light time equation (Equation G-1) 
is determined. The function p2(tp - At2) is expanded into a first-order 
Taylor series. The Taylor series for the function P(tj» - At2) at At2 = At2* 
is 


dPo(t - Ato) 

p2(tp - At2).= p2-(tj, - At2’) + 

d(At2) 


(At2 - At2') + . . . 


At2=At^ 


(G-5) 


where 


At2 = signal transmission time to be determined 
At2* = signal transmission time assumed 


dP2(tr> - At2) 
d(At2) 


P2 ■ V 
P2 

At2=At2' 


V = velocity vector of the vehicle at time t^, - At2*. 


Substituting the expression for P2 from Equation (G-1) into Equation (G-5) 
yields with some manipulation 


At. 


PsCtr - ^t2') ■ v(tj, - At2') 

Pp(t„ - Atp') + At; 

P2 (tr-At 2 ’) 


(G-6) 


e + 


P2(tp - At2*) * v(tj, - Atp') 

PpCtp - Atp') 


The scheme consists of iterating on Equation (G-6) until |Atp - Atp'l is 
sufficiently small: 
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Step 1: At 2 * =0. Obtain corresponding values of p 2 > n, v and p 2 - 


Step 2: At 2 = 


p2<tj,) 


P2(tj.) * V (tp) 


p2(tr,) 


Step 3: At 2 * = At 2 j value obtained from previous step. Obtain corres- 


ponding values of p 2 > r, v and-p 2 . 


Step 4: Compute At 2 using Equation (G-6). 

Step 5: Check if [At 2 ~ ^^ 2 '! 5 If it is sufficiently small, then 

the value- of t^ is tj, - At 2 - Final values of downlink range 
and range rate vectors plus their associated magnitudes are 
obtained from this corrected value of time. 

Step 6; If the value is unacceptable Steps 3> 4 and 5 are repeated 
until convergence is obtained. 


The uplink data is determined in a similar manner by iterating on the 
squation 


p^(t^ - At^') + 


At. 


Pl(ty - At-] ') • v(ty - At-] ' ) 
p^(t^ - At-]’) 


At-]’ 


c + 


P-](ty - At-]’) * V(ty.- At-]’) 


c + 


P-](t.y - At-]') ' vity - At-]')' 
Pi “ ^t-] ' ) 


(G7) 


ntil |At-] - At-] ’I is sufficiently small. Then the value for is t^ -Ati, 
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APPENDIX B.7 
VENT MODELS 


B. 7.1 GENERAL DESCRIPTION 

There are well over 50 purge, vent, and relief ports located throughout the 
Shuttle (e.g., fuel cell 0^, purge, avionics bay vents, H^O boilers, freon 
coolant loop flash evaporators, gas separator vent, APU dumps, ammonia dumps, 
urine dump vents, RCS fuel purge vents, etc.), as well as several sources of 
unscheduled leaking (e.g., cabin leakage, payload bay venting, etc.). Figure 
B-4 shows the location of major ports and vents. Of these sources the. APU and 
water spray boiler cause the largest perturbations to the Shuttle trajectory. 
During portions of ascent and deorbit, when the payload bay doors are closed, 
the topping flash evaporator and high' load flash evaporator steam vents are also 
major sources of trajectory perturbation. 

B . 7 . 2 ENVIRONMENTAL SIMULATION 

The following tables are from reference B.7-1. These tables list the vent source 
the location in Orbiter structural body coordinates (ref. B.7-2), the effluent, 
the magnitude of the thrust, the thrust direction in Orbiter body coordinates, 
and the duration of the vent. Table B.7-1 lists scheduled vents. Most vents 
listed vary with respect to time. Reference B.7-1 includes plots of venting 
magnitudes versus time. The APU venting thrust is also a function of pump shaft 
horsepower, while the hydraulic water spray boiler venting thrust is a function 
of pump shaft horsepower, water temperature, and time. Moreover, plume 
impingement effects are negligible for all but the three APU and three water 
spray boiler vents. A series of tables are given in reference B.7-1 that detail 
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TABLE SCHEDULED VENTS 


n 

I 

j\ 


Vent source 

Location , in 
_Io Xo_ - 


Effluent 

Thrust , 
Ibf 

X 

Vector 

y 2 

Duration 

Remarks 

ECLSS 











Cabin air (commode) 

602.0 

-105.5 

336.0 

Air 

Negligible 

0 

1 

0 

2 seconds, 1) to 7 
times per 24 hours 

Depends on number of 
people 

Water vapor from fecal 

material 

602.0 

-105.5 

336.0 

HgO 

Negligible 

0 

1 

0 

Same as above 

Same as above 

H- separator 

(FCP water) 

602.0 

-105.5 

336.0 

»2 

Negligible 

0 

1 

0 

Continuous 

, Flow rate 
2x10"'^ Ibm/sec 

Airlock depressur- 

ization 

Topping flash evapor- 
ator steam vent 

602.0 

-105.5 

336.0 

Air 

5.0 max. 

0 

1 

0 

5 minutes , three 
times per mission 

Dependent on fuel 
cell H 2 O production 

No. 1 

1505.7 

-128.lt 

301). 8 

H 2 O (steam) 

0.5 max. 

0 

1 

0 

Intermittent for 

rate. Supplements 

No. 2 

High load flash evap- 

1505.7 

-128.1) 

301). 8 

H 2 O (steam) 

0.5 max. 

0 

1 

0 

up to 10 hrs 

radiators onorbit. 
Used in conjunction 
with high load 
evaporators during 
ascent and deorbit. 

orators steam vent 

1390.3 

-113.8 

326.8 

H 2 O (steam) 

2.0 max. 

0 

1 

0 

Operational during 
ascent above 140 000 
ft and after pay- 
load bay door 
closure . 

Used in conjunction 
with topping evap- 
orator during ascent. 


Fuel Cell 
2 

H purge 

640.9 

105.4 

339.2 

H 2 (g) 

0.05 

0 -1 

0 

Normal mode - 2 
minutes per FCP 
every 8 hrs 

O 2 purge 

640.9 

- 105.4 

339.2 

O 2 (g) 

0.05 

0 1 

0 

Concurrent with H 2 
purge 


9 

S 

C5 

s 

1 E 

I 

> no 

" 

9 G5 
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TABLE Concluded 


Vent source 

Location, 
-J?o Y. 

in, 

^o_ 

Effluent 

Thrust, 

Ibf 

”x 

Vector 

- _L _I_ 

Duration 

Remarks 

Hydraulic water boiler 

vent 

No. 1 

1375.7 

16.3 

508.1 

H 2 O (steam) 

7 max. 

0 

0 -1 

6 minutes 

Postboost MPS dump; 
onorbit APU checkout. 
Ent.^’y preparation. 

No. 2 

1395.7 

16.3 

510.5 

H 2 O (steam) 

7 max.' 

0 

0 -1 

6 minutes 

No. 3 

1355.7 

16.4 

505.8 

H 2 O (steam) 

7 max. 

0 

0 -1 

6 minutes 


APU exhausts 

Left 1 

1327.9 

-15.9 

502.5 

N 2 (g) 

9.0 Ibf 

0 

0 -1 

Same as above 

Same as above 

Left 2 
Right 1 

1318.7 

1327.9 

-15.7 

15-9 

501.3 

502.5 

Ng (g) 
NH3 (g) 

average per 
vent , 30 Ibf 
peak total. 

0 

0 

0 -1 

0 -1 
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the force vectors and moment vectors of these vents. Tables B.7-II and B.7-III 
list contingency and failure vents, respectively. 

B.7.3 ERROR SOURCES 

Thrust uncertainties of the quantities listed in the tables are +5 percent 
for all vent forces except 

a. Plume impingement forces, which are +10 percent 

b. RGS and OMS pressure relief vent forces, which are +zero to -10 percent. 

These numbers are very optimistic and may be increased by a factor of 10 
(ref. B.7-3). 
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B.7.4 references 

B.7-1 Carminati, J. M.: Orbiter Vent/Thrust Reference Data For OFT Missions. 

MDTSCO Working Paper No. 1.4-7-204, January 4, 1978. 

B.7-2 Davis, L. P.: Coordinate Systems for the Space Shuttle Program. NASA 

TM X-58153, JSC-09084, October 1974. 

B.7-3 Os burn, R. K./FM8, Personal communication. 
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TABLE B.7-II.- CONTINGENCY VENTS 



Location, In 



Vent source 

Xo^ 


■ Zo_ 

Effluent 

ECLSS 
Urine dump 

620.0 

-105.4 

336.0 

Urine , 
sweat , 
biocide 

Potable water dump 

620.0 

-105.5 

342.4 

HgO ( 1 ) 

Fuel Cell 
Water relief 

631.1 

105.4 , 

339.2 

HgO ( 1 ) 

Reactant regulator 
Vent relief 

6H0.9 

640.9 

-105.4 

-105.4 

339.2 

339.2 

O 2 (g) 
H 2 (g) 

Hydraulics water boiler 

vent 

No. 1 

No. 2 

No. 3 

1375.7 

1395.7 

1355.7 

16.3 

16.3 

16.4 

508.1 

510.5 

505.8 

HgO (g) 
HpO tg) 
HgO (g) 

APU exhausts 
Left 1 
Left 2 
Left 3 

1327.9 

1312.7 

1327.9 

-15.9 

-15.7 

15.9 

502.5 

501.3 

502.5 

N 2 (g) 
H 2 (g) 
NH 3 (g) 


Thrust, 
Ibf 

X 

Vector 
Y 2 

Duration 

Remarks 

1 .275 

0 

1 

0 

4.5 hr/week-mission 
> 7 days 

No dump for 7-day 
mission even with 
double failure 

1.275 

0 

1 

0 

1 .5 hr/day 

Used if failure 
in flash evaporator 
or for thermal under- 
load 


< 0.01 

. 0 -1 

0 

Can 

be 

continuous 

Occurs only if 
ECLSS cannot accept 
PCP water 

0.05 

0 1 

0 

Can 

be 

continuous 


0.05 

0 -1 

0 






7 max. 

0 

0 -1 

6 minutes 

Checkout of 

7 max. 

0 

0 -1 


hydraulic system 

7 max. 

0 

0 -1 




9.0 Ibf 

0 

0 -1 

6 minutes 

Checkout of 

average per 

0 

0 -1 


hydraulic system 

vent 30.0 

0 

0 -1 




Ibf peak 
total. 


r" 

' ^ 

■ rn 
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Looation, in . 


Vent source Effluent 


RCS pressure relief 





He+ 

Fwd oxidizer 

285.5 

40.6 

361.4 

NpOi, 

MMH 

Fwd fuel 

285.5 

-40.6 

361 .4 

Aft oxidizer (RHS) 

1507.7 

128.1 

510.1 

NpOh 

MMH 

Aft fuel (RHS) 

1509.11 

65.7 

529.9 

Aft oxidizer (LHS) 

1507.7 

-128.1 

510.1 

NpOi, 

Aft fuel (LHS) 

1509.1* 

-65.7 

529.9 

MMH 

OMS pressure relief 
Oxidizer (RHS) 

1507.7 

128,1 

510.1 

He+ 

NjOi, 

Fuel (RHS) 

1509.11 

65.7 

529.9 

MMH 

Oxidizer (LHS) 

1507.7 

-128.1 

510.1 

NpOi, 

MMH 

Fuel (LHS) 

1509.4 

65.7 

529.9 


OMS pressure relief 
(PLB) 

Oxidizer 

Fuel 

1295 

1295 

105 

-105 

331.2 

331.2 

He+ 

NgOi, 

MMH 

APU fuel pump seal vents 
APU 1 
APU 2 
APU 3 

1343.3 

1349.9 

1343.3 

-107.8 

-108.3 

107.8 

413.3 

413.6 

413.3 

NpHi, 

NgHi, 

NgHj, 

ECLSS 
O 2 vents 

602.0 

-105.5 

336.0 

Og (g) 

N 2 vent 

602.0 

-105.5 

336.0 

Ng (g) 
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TABLE B.7-III.- Concluded 


Vent source 

Location, in 
Yo 

_Zo_ 

Effluent 

Thrust, 

Ibf 

~X~ 

Vector 

y z 

Duration 

Remarks . 


Cabin overpressure 

571.0 

-67.7 

394.4 

Air 

53.0^’® 

-1 

0 

0 

Flow stops at crew 

Opens at 15.5 psia 


Relief valves'^ 

Power reactant storage 
and distribution 

576.0 

-?lt.2 

394.4 

Air 

53.0^'® 

-1 

0 

0 

initiated shutdown. 

cabin pressure. • 


H2 relief vent 

622.2 

105. 4 

342.5 

H2 (g) 

60^ 

130® 

0 

0 

-1 

-1 

0 

0 

Total venting time 
(each tank) is 17 
hrs; 20 min. in 
cycles < 10 sec 

Occurs due to tank 
overheating; heaters 
■failed on most likely 
cause. 

Hg multiple-tank relief 

vents 

990.0 

105.11 

352.2 

H2 (g) 

1071° 

0 

-1 

0 

Same as above 

Same as above 


Og relief vent 

760.0 

-105.11 

355.5 

■O2 (g) 

17b 

0 

1 

0 

Total venting time 
(each tank) in I6 
hrs; 2.4 hrs in 
cycles < 100 sec 

Occurs due to tank 
overheating; heaters 
failed on most likel, 
cause. 

V 

Og 12 tank relief 
Fuel cell 

760.0 

-105.11 

355.5 

O2 (g) 

600^ 

0 

1 

0 

Same as above 

Same as above 

0 0 

H2 purge 

6110.9 

105.4 

339.2 

H2 (g) 

0.3 

0 

-1 

0 

Flow stops at crew 
initiated isolation 

Failure of FCP 
reactant regulator 

■n ^ 
tj fi5 
Q 2 

Og purge 

Hydraulic water boiler 
vent 

6110,9 

-105.4 

339.2 

O2 (g) 

0.5 

0 

1 

0 

Same as above 

Same as above 

0 5 ? 
TO F* 

lO 

c :©* 

No. 1 

1375.7 

16,3 

508.1 

HgO (g) 

7 max. 

0 

0 

-1 

' No venting if 


> Q> 
r" TO 
3 

No. 2 

1395.7 

16.3 

510.5 

HgO (g) 

7 max. 

0 

0 

-1 

hydraulic system 


No. 3 

1355.7 

16.4 

505.8 

H2O (g) 

7 max. 

0 

0 

-1 

fails 


^ m 


^Plume impingement within PLB will change magnitude and direction of thrust vector. 

“Two tanks relieving simultaneously. 

Three tanks relieving simultaneously through the forward vent (Xq 822.2) and nine relieving simultaneously through the multiple tank 
relief vents, 30 day mission only. 

“Twelve tanks relieving simultaneously, 30 day mission only. 

^Regulator on line failures resulting in flow through both valves. 
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TABLE B. 7 -IV.- RCS AND QMS PRESSURE RELIEF VENTS 




Location, 

in. 

Direction cosines 


RCS 



^o_- 

i • 

J 

k 

Fwd 

oxidizer 

285.5 

40.6 

361.4 

- 0.4520 

0.8103 

0.3729 

Fwd 

fuel 

285.5 

- 40.6 

361 .4 

- 0.4520 

- 0.8103 

0.3729 

Aft 

oxidizer (RHS) 

1507.7 

128.3 

510.1 

0.9986 

0.0379 

0.0366' 

Aft 

fuel (RHS) 

1509.4 

65.7 

529.9 

0.9986 

0.0379 

0.0366 

Aft 

oxidizer (LHS) 

1507.2 

-128.3 

510.1 

0.9986 

-0.0379 

0.0366 

Aft 

fuel (LHS) 

1509.4 

- 65.7 

529.9 

0.9986 

-0.0379 

0.0366 



Location, in. 

Direction cosines 

OMS 

„Jo_ 


2 ^ 

i 

.1 

k 

Oxidizer (RHS) 

1507.7 

128.0 

510.1 

0.9986 ' 

. 0.0379 

0.0366 

Fuel (RHS) 

1509.4 

65.7 

529.9 

0.9986 

0.0379 

0.0366 

Oxidizer (LHS) 

1507.2 

-128.1 

510.1 

0.9986 

-0.0379 

0.0366 

Fuel (LHS) 

1509.4 

- 65.7 

529.9 

0.9986 

-0.0379 

0.0366 
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APPENDIX B.8 
VEHICLE BENDING 


B.8.1 GENERAL DESCRIPTION 

In most simulations the Space Shuttle is treated as a rigid body. In analyzing 
the effects of external navigation aids on the navigated state, this is a good 
approximation. However, when analyses of guidance, of flight control, or of the 
effects of internal navigation aids, such as the IMU's or RGA’s, are desired 
vibration and bending modes of the vehicle are relevant for realistic modeling. 
This appendix gives a brief overview of the theory of elastic motion as employed 
by the space vehicle dynamics simulation (SVDS) . For more detail see references 
B.8-1. and. B.8-2. 


It is assumed that the equations of motion of a system can be derived from a 
Lagrangian function defined by the kinetic energy function T, the virtual work 
of the external forces 6W, , the strain energy function U, and the Rayleigh 
dissipation function R. These can be written as 


T = l/2(qj^ qgT qg'^] 


1 


<1r 

%R %B ^BS, 


<3B 

^SR ^SB Mgs 


9S 


T T 

= <Sqg Qg + 6qg Qg 

D = (qg^ Fgqg + qg^ Pgqg)/2 
R = (Qb'^^B^B 

where qg, qg, and qg are respective sets of pseudorigid body generalized 
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coordinates, elastic body generalized coordinates defined by a set- of base 
modes, and linear springs mass generalized coordinates. is the modal 

coupling mass matrix, is the generalized force acting on the body, is 

the modal stiffness matrix, and Rj^. is the Rayleigh dissipation function. The 
subscripts refer to pseudorigid body coordinates (R), elastic body coordinates 
(B), and linear spring coordinates (S). In terms of the above quantities 
Lagrange's equations can be written as 

d 3T 3U 3R ^ 

dt 3q^ 3qj^ 3q^ 3q^^ ^ 

subject to the constraint 


%R ^R %B ^B %S 'is = 

The following set of equations result in 

^R'IR ^B'iB + “RS'iS = ° 

.. .. .. ^ 

^R'iR ^B'iB ^BS*iS ^B*iB * ^B'iB = Qb *■ ^BR^RR % 
Msrqr + Msbqb + Mgsds + Rs'is + ^s'is = -^srMrr ^Qr 
These equations can be rewritten as 

®B = ®BB % -^BR^RR ^ ®R '^B'iB " ^B^^B 
+ %S®SS ^ ^%R%R ^ Qr 

•^S ■ "^SS ^^SI^RR ®R * ^S'is 
-1 

•iR = -%R (MpB'JB + ^S'ls) 
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where Gsb = Msg - MsRMgg-l Mpg 
T 

^BS = ^SB 

®ss = %s “ %r”rr ^ %S 

<^BB = Mbb - MbrMbr ^ Mrb - GgsGss ^ Gsg 

The constraint equation yields the following two equations that, with equations 
(2) through (4), completely define the pseudorigid body motion, 

-1 

9r = “Mrr * ^RS9S^ ^5) 

'^R “ "^RR ^ 

The actual displacement of the point of interest is given by 


P = <t>Q - <t>RqR + 


Where <j) is the matrix of base mode shapes and can be broken down into two 
subraatrices, 4 >r> of rigid body base mode shapes, and 4>g, elastic body base 
mode shapes which include both elastic body and linear spring modes. The 
matrix, <i), has the properties that at some time, T, 




%R 

0 


0 

I 





and 


T 

<j> B4> 



0 


T T 


M(t) 


F(T) 


R(t) 
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Where A(t) is the panel point inertial matrix; K, the panel, point stiffness 
matrix; B, the panel 'point damping matrix. The notation, C denotes a 
diagonal matrix. oj.j. are the natural frequencies corresponding to elastic base 
modes; is the critical damping ratio for the elastic base modes. The 

velocity and acceleration of the panel point are given by 

P = <j>q 
and p = 4>q 

Reference B.8-3 gives data for various vehicle configurations and various 
missions for use in studying vehicle bending. 

Although the navigation principal function runs at a low rate of 0.25 Hz, there 
are other navigation related functions that run at higher frequencies; e.g., the 
user parameter processor at 6.25 Hz and the attitude processor at 12.5 Hz. The 
spectrum of vehicle bending frequencies starts at 2 Hz and continues to 
infinity. The higher rate functions, which employ IMU and RGA data, will be 
affected in a systematic way, irtiereas the lower rate functions, navigation and 
redundancy management, will be affected in a more or less random manner. 
Furthermore , modeling errors will occur due to the finite cutoff of the bending 
spectrum when using a digital simulation of this essentially analog process. 
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B.8.2 REFERENCES 

B.8-1 SVDS Program Subroutine Library, Vol. I, Subroutines A-C. JSC IN 73- 
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B.8-2 Flanders, H. A.; Flexible Body Capability in the SVDS Simulation. . 
TRW I.C. 6531.2-66, October 31, 1972. 

B.8-3 Space Shuttle Flight Control System Data Book, Vol. I, Integrated 
Vehicle. Rockwell International SD73-SH-0097-1E, April 1976. 
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APPENDIX B.9 
CLOCK AND TIMING ERRORS 

L9.I GENERAL DESCRIPTION 

The Space Shuttle master timing unit (MTU) provides time and synchronization 
frequency to various users. Greenwich mean time (GMT) and mission elapsed time 
(MET) are provided to all computers, displays, recorders, and telemetry (PCM). 
The GMT that is provided is universal coordinated time (UCT) (see appendix B.10 
for a discussion of UTC) . 

Synchronization frequencies are^ provided, to one-way Doppler systems, PCM master 
units, static inverters, and the time display (event -timer). , , 

Power is supplied redundantly to the MTU through two separate circuit breakers 
from two separate essential power busses. A switch is provided to manually 
select one of two redundant oscillators or to enable the MTU to automatically 
select the oscillator based on internal BITE (built-in test equipment) 
information. 

B.9.2 ERROR SOURCES 

The master timing unit performance is as follows: 

B . 9 . 2 . 1 Master Timing Unit Performance (ref. B.9-1) 

a. Frequency accuracy (offset): +7x10“® 

(can be adjusted to within +1x10”^) 

b. Frequency stability (drift): +1x10“^ per day 

c. Time accuracy: +10 milliseconds bias 

one-eighth millisecond quantization 
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Mission rules permit oscillator frequency offset to be as much as +7x10~® 
at the start of any 7-day mission. Greater offset will require that the 
oscillator be adjusted before flight. Periodic adjustments at the launch site 
are anticipated without removal of the MTU- 

d . Time accuracy - GMT will be- set before launch to within +1 millisecond 

referenced to U.S. Naval Observatory time. GMT can then drift off as much 
as +10 milliseconds at any time during the mission. A frequency offset of 
+8x10“*^ will result in a time error of approximately +10 milliseconds in 24 
hours. It is anticipated that a GMT update every 24 hours may be required to 
keep the MTU within the 10 millisecond limit. GMT initialization is accomp- 
lished by the ground computer launch processing system (LPS) via the 
umbilical and onboard computer. In flight the time delta is determined by 
the MCC ground computers and the GMT' is then corrected via the uplink. It 
can also be updated by voice to the crew who can then input it by the 
keyboard through the computer to the MTU. 
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B.9.3 REFERENCES 

B.9-1 Procurement Specification for Master Timing Unit. Rockwell International 
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APPENDIX B.10 
TIME 


B.10.1 GENERAL DESCRIPTION 
B. 10.1.1 Ephemeris Time 

Ephemeris time is the fundamental unit of time and is based on the period of 
orbital motion of the Earth about the Sun; thus, the time unit is in- years. An 
ephemeris second is a subdivision of a period, or, in exact terms, it is 31 556 
925.97^7~^ of the tropical year for 1900. A tropical year is the time required 
for the Earth to rotate one revolution about the Sun passing through the point 
of vernal equinox. 

B.10. 1.2 Atomic Time 

Atomic time is based on the time required to accumulate a fixed number of 
vibrational periods of transition frequency of atoms under a prescribed 
condition. The frequency of the atomic transition and, thus, the period of the 
frequency or any multiple thereof is so constant that by comparision with the 
ephemeris time it takes at least a century to be certain of any atomic time 
variation . 

The Twelfth General Conference of Weights and Measures, which met in October 
1964, designated "an atomic or molecular frequency standard to be employed 
temporarily for the physical measurement of time" and further declared "that the 
standard to be employed is the transition between the hyper fine levels F = f,mf> 

= 0 and F = 3, mp = 0 of the fundamental state 2S-j/2 atom of Cesium 133, 

not perturbed by external fields and that the value 9 192 631 770 Hertz is 
assigned to the frequency of this transition." F, ra^, and 23^/2 represent 
numbers that identify the quantum states and energy levels associated with a 
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particular atomic transition and hence, vath a particular spectral line (or 
frequency) of emission. 

B. 10. 1.3 Universal Time or Greenwich Mean Time 

Universal time (UT) is based on the Earth's period of rotation about its axis; 
thus, the fundamental time unit is in solar days. Since apparent solar days vary 
in length with the season, mean solar time was devised, and the tabulation of 
the difference between the apparent value and the mean value is called the 
equation of time. Because of its definition, UT is dependent upon the obser- 
ver’s position; thus, a reference based on the prime meridian at Greenwich, 
England was established. UT is therefore synonymous with GMT. The following 
are important points concerning UT (ref. B.10-1). 

a. Different versions of universal time - The value of UT at a particular 
instant must be determined from astronomical observations. This observed 
universal time (UTO) is calculated by the USNO and is based on observations 
of star transits across the meridian“lihe passing through Washington, D.C., 
and Richmond, Florida. This version of UT is subject to a number of irregu- 
larities caused by variations in the rotation rate of the Earth. A more 
uniform UT (UT2) is obtained by correcting UTO for the seasonal variation 
variation variation in the rotation rate of the Earth and for polar drift. 

The correction formula is 

UT2 = UTO + S + P 


where 


S is correction for seasonal variation 
P is correction for polar drift 
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The precise value of UT2 cannot be known until after the astronomical data 
have been processed. The coordinated universal time (UTC) system provides 
a measure of universal time that is not subject to this delay. It is a 
convenient, up-to-date approximation of UT2 that is based on an atomic clock 
constrained to be in close agreement with UT2. 

b. Time system maintained by USNO and NBS - The National Bureau of Standards 
(NBS), through its radio station, WWV, transmits time signals monitored by 
its UTC clock. The stations of the National Aeronautics and Space 
Administration (NASA) spaceflight tracking and data network (STDN) use the 
WWV signals to time tag the observed- Doppler and range data, and the clocks 
onboard spacecraft are calibrated prelaunch with WWV. In addition, all 
time- dependent calculations, both onboard and on the ground, use UTC as the 
basic dependent variable. 

Under the old UTC system the rate of the atomic clock was changed by 
offsetting its output frequency, which is based on the natural resonant 
frequency of Cesium atoms, by an amount that made the rate of the clock 
nearly identical with the average rate of rotation of the Earth with all 
known seasonal and polar variations- removed (UT2). The amount of rate 
difference was the offset of the UTC with respect to atomic time. On 
occasion, whenever the adopted frequency offset would clearly allow the UTC 
to diverge from UT2 by more than 0.1 second, the UTC clock of NBS would make 
instantaneous step adjustments. These periodic time corrections are known 
as step- time adjustments. The offset and step- time adjustments were made at 
the same time and by the same amount throughout the world by international 
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coordination through the Bureau International de 1‘Heure (BIH); thus, the 
term coordinated universal time. 

c. Problems with the old system - The variable offset method used in the pre- 
1972 OTC system prevented frequent step- time adjustments and answered 
astronomical and navigation needs each time an offset change was introduced; 
however, the length of the broadcast time interval (the UTC second) changed. 
Although the pre-1972 system satisfied those vdio required time based on the 
rotation of the Earth, other users needed time intervals that did not change 
and that were more directly related to atomic standards. 

In 1968 CCIR (International Radio Consultant Committee) recognized the need 
to make the broadcasted second conform to the international definition of 
the second (the duration of 9 192 631 770 periods of radiation corresponding 
to the transition between the hyperfine levels of the ground state of the 
Cesium 133 atom). A working group was established to study the problem and 
its recommendations were adopted in February 1971 and implemented in January 
1972. 

d. The improved UTC system - The improved UTC system eliminates the changing 
frequency offset and increases the size of the step-time adjustments from 
0.1 to 1 second, which is called a leap second. Therefore, the UTC clocks 
are driven at the atomic clock rate but, by making periodic leap second 
adjustments, are synchronous with UT to a tolerance of +0.7 second. Figure 
B-5 shows the difference between the old UTC system prior to 1968 and 

the new UTC system in use at that time. Station WWV changed its clock 
time from the old UTC system to the improved UTC system by a time step 
that occurred on December 31 » 1971, at 23 hours 59 minutes 60.107600 
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seconds (old UTC) . At that instant the date became January 1, 1972, 0 
hours 0 minutes 0 seconds exactly (new UTC). To users such as navigators, 
surveyors , and geodesists ^o need UT1 (UT corrected for polar variation 
only), the time difference between UT1 and UTC is transmitted in code to a 
precision of 0..1 second. 



Figure B-5.- Recent behavior of TAI-UTC. 
e. Ephemeris time and its relation to UTC - The measure of ET is obtained by 
coraparision of observed positions of the Sun, Moon, and planets with their 
corresponding ephemerides . The ephemerides of the Sun and Moon are needed 
for perturbation computations, targeting, and alinements and are maintained 
in ET. Since timekeeping, both onboard and on the ground, is in UTC the 
difference (At) between ET and UTC must be available. Since atomic clocks 
vary at the same constant rate as ephemeris time, and the relationship 
between atomic time and UTC is known, the At is easily computed. As of 
July 1977 the following formulae apply: 
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Hence , 


ET = TAX + 32.18 seconds (1) 

TAX - UTC = 16.0 seconds (2) 

At = ET - UTC = 48.18 seconds 

where TAX is XnternatlonaX atomic time and UTC is NationaX Bureau of Standards 
UTC. Equation (1) is constant. 

B.10.2 Error Sources 

The clocks associated with the ground tracking equipment are synchronized by 
time signals received via the LORAN-C networks. These networks are in turn 
synchronized to a master UTC clock of the United States Naval Observatory 
(USNO-MC) to an accuracy of 5 ys. Some LORAN-C networks, such as the U.S. East 
Coast, North Atlantic, and Norwegian Sea chains are synchronized to the USNO-MC 
within 0.5 ys, but the Central and Northwest Pacific may be off as much as 5 ys 
of the USNO-MC. 

The S 3 mchronization of ground station clocks to each other was calibrated during 
Apollo 16 by the use of interstation time synchronization tests. The results 
showed that the largest dispersion of any reference receiver combination was 5 
ys with an average dispersion of less than 2 ys. 

The use of portable clocks, and the complete implementation of the Defense 
Satellite Communication System (DSCS), known as Comm-Sat, will improve the 
synchronization of the LORAN-C network to the USNO-MC and the site-to-site 
synchronization to less than 1 ys worldwide, but these systems are currently not 
fully operational and would require some changes in ground station 
synchronization procedures. 
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APPENDIX B.11 
ASTRODYNAMIC CONSTANTS 

This section provides the basic astrodynamic constants for use in error analysis 
and mission simulations related to Space Shuttle development. The data included 
here are provided to facilitate the use of consistent numerical values 
throughout the Johnson Space Center (JSC). 

The values of the constants presented here are taken primarily from reference 
B. 11-1 . 

B . 1 1 . 1 CONVERSION FACTORS 

The conversion factors that are used to convert between systems of units are 
presented in the table B.11-I. These are the minimvim number of factors required 
for conversion between the systems of units commonly in use at JSC. 

TABLE B.11-I.- CONVERSION FACTORS 


Original unit 

Earth radius 
Foot 

Nautical mile 
Pound- mass (Ibm) 
Hour (time) 
Radian 


Desired unit 
Meter 
Meter 
Meter 
Kilogram 
Second (time) 
Degree 
Arc-second 


C onversion factor 

6.378160 X 10 

-1 

3.0il8 X 10 
1.852 X 10^ 
4.5359237 x 10"^ 

3.600 X 10^ 

180/tt 

3.600 X 10^ 




Degree 
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b.11.2 constants and parameters 

B . 1 1 . 2 . 1 Earth Constants 

B.11.2. 1.1 Equatorial Earth radius (Gravitational) 

1 e.r. = 0.6378160000000000+007 + 0.5000000000000000+001 m (ref. B.11-1) 

1 e.r. = .2092572178477690+008 + .1640419947506562+002 int. ft 

1 e.r. = .1000000000000000+001 + .7839251445557967-006 E.r. 

1 e.r. = .3443930885529158+004 + .2699784017278618-002 n. mi. 

B.11.2. 1.2 Flattening 

1/f = 0.2983500000000000+003 + 0.4000000000000000-001 (ref. B.11-1). 

f = 0.3352891869237217-002 

B.11.2. 1.3 Gravitational Parameter (GMg = lig = li Earth): 

Vig = 0.3986012000000000+015 + 0.4000000000000000+009 m^/sec^ (ref. B.11-1) 

Vig = .1407646853278542+017 + .1412586668859544+011 (int. ft)^/sec^ 

3 2 

Ug = .1990936344044169+002+ .1997923081059631-004 (E.r.) /hr 

yg = .6275027808522208+005 + .6297048587432459-001 (n. mi.)^/seo^ 

B.11.2. 1.4 Mass of the Earth 
Mg = 0.5973343323842350+025 kg 

Mg = 0.1316896781981220+026 Ibm 

B.11.2. 1.5 Angular Velocity of the Earth's Rotation With Respect to the Vernal 
Equinox 

By differentiating Newcomb’s expression of the right ascension of the fictitious 
mean Sun (ref. B.11-2), one obtains 

% = 9856473354 + 2.1203 x 10"® Ty deg/d (1) 
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where Ty is Julian centuries of 36 525 days counted from 1900 January 0, 

12h UT and Ry is the time rate of the right ascension. 

Adding to this the. time rate of the hour angle, 360 deg/d 

-8 

cop = 360.9856473354 + 2.1203 X 10 Tu deg/d 
or 

cop = 7.292211585458 X 10"^ + 4:283 X 10 Tu rad/sec (2) 

gives the rotation rate of the Earth with respect to a processing vernal 
equinox. Tabular data for each year to the year 2000 is given in table B. 11-11. 

The motion defined in. equation (2) is the result of the spin of the Earth and 
motion of the vernal equinox (precession) . Because the latter motion takes 
place in the ecliptic, the equitorial component, 'or the general precession in 
the right ascension, m, must be used to compute the inertial (nonprecessing) 
Earth rate. 

Wj = o)p “ m . 

From reference B.11-3 

-8 -12 2 

, m = 0.00003505464 + 2.1253 x ,10 Ty + .92 x 10 Ty deg/d 
The angular velocity of the Earth's rotation with respect to an inertial 
coordinate system is therefore 

-8 

= 360.9856122808 - .0050 x 10 Ty deg/d 
or 

coj = 7.29211514646 x 10"^ - .010 x 10“^^ Ty rad/sec 
This rate is constant to 13 significant digits for the rest of this century. 
Tabular values are given by year until the year 2000 in table B. 11-11. 
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TABLE 

A 


Year 

1960. 

1961 . 
1962. 

1963. 

1964. 

1965. 

1966 . 

1967. 

1968. 

1969. 

1970. 

1971. 

1972. 

1973. 

1974. 

1975. 

1976. 

1977. 

1978. 

1979. 

1980. 
1981. 
1982. 

1983. 

1984. 

1985. 

1986. 

1987. 

1988. 

1989. 

1990. 

1991 . 

1992. 

1993. 

1994. 

1995. 

1996. 

1997. 

1998. 

1999. 

2000. 


OF POOR QUALITY 

79FM5 

B. 11-11.- THE EARTH'S ANGULAR ROTATIONAL VELOCITY WITH RESPECT TO 
PRECESSING EQUINOX (Wp) AND INERTIAL EQUINOX (uj) VERSUS YEAR 


0 . 72921 1 5854836980-004 
. 7292 1 1 585484 1 263-004 
.72921 15854845546-004 
.72921 15854849829-004 
.7292115854854112-004 

.7292115854858395-004 

.7292115854862678-004 
.7292115854866961-004 
.72921 15854871244-004 
.72921 15854875527-004 
.729211 58548798 1 0-004 
.7292115854884093-004 
.7292115854888376-004 
.7292115854892659-004 
.7292115854896942-004 

.72921 15854901225-004 
.7292115854905508-004 
.7292115854909791-004 
.729211 58549 1 4074-004 
.72921 15854918357-004 
.72921 15854922640-004 
.7292115854926923-004 
.7292115854931206-004 

.7292115854935489-004 
.7292115854939772-004 
.72921 15854944055-004 
.7292115854948338-004 
.7292115854952621-004 
.7292115854956904-004 

. 7292 1 1 585496 1 1 87-004 
.7292115854965470-004 

.7292115854969753-004 

.7292115854974036-004 
. 7292 1 1 58549783 1 9-004 
.729211 5854982602-004 

.7292115854986885-004 

.7292115854991168-004 

.7292115854995451-004 

.7292115854999734-004 
.72921 15855004017-004 
.7292115855008300-004 


Uj 

0.72921151 46459400-004 
.7292115146459390-004 
.7292115146459380-004 
.72921 15146459370-004 
.7292115146459360-004 
.7292115146459350-004 
.7292115146459340-004 
.7292115146459330-004 
.72921 15146459320-004 
.7292115146459310-004 
.7292115146459300-004 
.7292115146459290-004 
.72921151 46459280-004 
.7292115146459270-004 
.7292115146459260-004 
.7292115146459250-004 
.7292115146459240-004 
.7292115146459230-004 
.72921 15146459220-004 
.72921 15146459210-004 
.72921 15146459200-004 
.7292115146459190-004 
.72921151 46459 1 80-004 
.7292115146459170-004 
.7292115146459160-004 
.7292115146459150-004 
.7292115146459140-004 
.7292115146459130-004 
.7292115146459120-004 
.7292115146459110-004 
.7292115146459100-004 
.7292115146459090-004 
.7292115146459080-004 
.7292115146459070-004 
.7292115146459060-004 
.7292115146459050-004 
. 72921 1 5 1 46459040-004 
.7292115146459030-004 
.72921151 46459020-004 
.7292115146459010-004 
.7292115146459000-004 
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^AGE fs 

QUALfTY 

B.11.2.2 Fischer Earth Model 
The following constants describe the Fischer Earth model (I960), which is used 
for the location of radar stations for the Space Shuttle. 


OF POOR 


B.11.2.2.1 Equatorial Earth radius 
a = 0.6378166000000000+007 m 
a = .20925741^16981627+008 int. ft 

a = . 1000000940710173+001 E.r. 

a = .3443934125269978+04 n. mi. 

B.11.2.2.2 Flattening 

f = flattening = 1 - b/a 
f = 1/298.30 = 0.3352329869259135 X 10 


B.11.2.2.3 Polar Earth radius 

b = 0.635678428360717+007 m 
b = .2085559148165061+008 int. ft 

b = .9966486076873435+000 E.r. 

b = .3432388922034075+004 n. mi. 

B.11.2.2.4 Eccentricity of Ellipsoid 
e = (a^ - b^)/a^ 
e = 0.8181333401693114-001 
e^= 2f - f2 

e = 0.6693421622965943-002 
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B . 1 1 . 2 . 3 B asic Equivalents and Conversion Factors 

1 int. ft = 0.3048000000000000+000 (exact) m (ref. B.11~1) 

1 n. mi. = .1852000000000000+001 (exact) km (ref. B.11-1) 

1 km = .5399568034557235+000 n. mi. 

1 m = .3280839895013123+001 int. ft 

1 n? = .1076391041670972+002 int. ft^ 

1 n. mi. = .6076115485564304+004 int. ft 

1 rad = .5729577951308233+002 deg 

1 deg = .1745329251994329-001 rad 

1 Ibm = .4535923700000000+000 kg (exact) (ref. B.11-1) 

1 kg = ..2204622621848776+001 Ibm 

1 km = .1567850289111593-003 E.r, 

1 E.r. = .6378160000000000+004 km (exact for scaling) 

1 U. S. Survey ft = 0.3048006096012192 m 
1 int. stat. mi. = 5280 (exact) ft 
1 Ibf = 32.174048556 Ibm int. ft/sec^ 

TT = 0.3141592653589793+001 (ref. B.11-1) 

1 E.r. /hr = 0.5812700495771362+004 int. ft/sec 
1 int. ft/sec = 0. I720370765236369-003 E.r. /hr 
1 m/sec = 0.5644261040801736-003 E.r. /hr 
1 E.r. /hr = 0.1771711111111111+004 m/sec 

S.11.2.4 General Constants 
1.11.2.4.1 Astronomical Unit 

AU = 0.1495978930000000+012 + 0.5000000000000000+004 m (ref. B.11-1) 
AU = .4908067355643045+012 + .1640419947506562+005 int. ft 
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AU = .23^5470997905352+005 + -7839251445557967-003 S^.r. ■ 

AU = .8077640010799136+008 + .2699784017278618+001 n, mi. 

B. 11. 2. 4. 2 Velocity of Li^t in a Vacuum: 

c = 0.2997925000000000+009 + 0.300000000000000+003 m/sec (ref. B.11-1) 
c = .9835711942257218+009 -+ .9842519685039370+003 int. ft/sec 

G= .1692107128074554+006+ .1693278312240521+000 E.r. /hr 

c = .1618750000000000+006 + .1619870410367171+000 n. mi. /sec 

B. 11.2.4.3 Gravitational Parameters for the Sun: 

^ Sun = 0.1327124990000000+021 +0.1500000000000000+014 m^/sec^ (ref. B.11-1 
^ Sun = -4686697671960888+022 + .5297200008223289+015 (int. ft)^/sec^ 

U Sun = -6628734122426762+007+ -7492211553969865+000 (E.r. )^/hr^ 

Sun = -2089242635906454+011 + .2361393220287172+004 (n. mi.)^/sec^ 
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APPENDIX C 

SUBROUTINE LISTINGS AND DESCRIPTIONS 
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APPENDIX C, 1 
RANDOM NUMBER GENERATOR 

The subroutines to be used for generation of random numbers in simulation 
programs are described in this appendix. 

Two characteristics are desired. The numbers should be random in the sense that 
no patterns exist that affect the simulation, but they should follow a deter- 
ministic sequence in order to facilitate debugging any program modifications. 
These subroutines have the desired characteristics . 

The subroutines described here use the congruential method, which is discussed 
in reference C-1 . 

The specific implementation described here is for the Uni vac 1108, but may be 
readily adapted for any similar machine. 

The first step in the generation of the rando 

Wf 

to generate a random number from a set uniformly distributed between zero and 

35 

one. Let N be an (odd) integer, less than 2 -1; an initial value of unity is 

often used. The new value of N is given by the least significant 35 bits of the 

15 

product of the old value and 5 . The least significant 27 bits of N are then 

combined with an exponent and converted to a normalized floating point number, 

A. The output "A** is uniformly distributed over the interval zero to one. 

The next step is conversion of the uniformly distributed random numbers to 
normally distributed random numbers. This transformation is accomplished by 
using two uniformly distributed numbers, A^ and in the following equations 


^.numbers is to use subroutine RDM 
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= (-2 In cos (21T A^) 

1/2 

= (-2 In sin (2ir A^) 

These equations are implemented in subroutine RN2S. Every other call to RN2S 
results in two calls to RIM to generate A^^ and A^ and the corresponding 
and X^; then is returned. On alternate calls X^ is returned. The 

particular version of RN2S, listed on page C-8, has several interesting 
features- By setting RNUM and X2 equal to zero the random numbers generated by 
RN2S are all zero; thus, a particular error model can be made perfect for 
checkout purposes. Secondly, use of the limit check value, SIGMA, will ensure 
that nominal systems, i.e., errors limited to Id or 3 ^, will not be multiplied t 
large random numbers. If a number larger than SIGMA is chosen, it is thrown 
away and a new number is generated. In this manner the shape of the Gaussian 
curve is preserved; i.e,, no spikes at the end points. 
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CD * * * '•< A * A Yr‘ * * * A ;V A * lA * !! A A A A A A A A A A A .. A -• 

CD 

CD PROGRAMMER AND DATE 
CD . T. P. KINNEY 

CD DOCUMENTATION 
CD - , J . R. BRi)UN 

CD ■ ‘OCTOBER, 1969 

CD 

CD PURPOSE . • 

CD TO GENERATE SAMPLES FROM A UNIFORM DISTRIBUTION ON T'C 

CD INTERVAL FROM ZERO TO ONE. 

CD 

CD USAGE 

CD CALL RDMCI, A) 

CD ' 

CD DESCRIPTION OF PARAMETERS 
CD 

CD INPUT 

CO CALLING SEQUENCE 

CD ■ I •> INTEGER STARTER IS INPUT 

CD COMMON 

CO ■ NONE 

CD CARD 

CO NONE 

CD TAPE 

CD NONE 

CD 

CD OUTPUT' 

CD , CALLING SEQUENCE 

CD I - RESULT IS OUTPUT 

CD I A - NORMALIZED FLOATING POINT SINGLE PPECISION NUMJER, 

CD ' UNIFORMLY DISTRIBUTED ON INTERVAL FRvM ZERO TO ONl 

CD COMMON 

CD NONE 

CD CARO 

CD ■ , NONE 

CD PRINT 

CD NONE 

CD TAPE 

CD none 

CD 

CD REMARKS AND RESTRICTIONS 

CD . I .NE. 0 AND NOTE THAT VALUES DO REPEAT AFFEP 

CD NUMBERS HAVE BEEN GENERATED AND THAT ANY GENERATED VALUE 

CD IS A FUNCTION OF THE NUMBERS USED TO INITIALIZE THE 

CD ■ PROCEDURE 

CD 

CD SUBROUTINES REQUIRED 
CD NONE 

CO 

CD METHOD 

CD THE FIXED POINT NUMBER I IS MULTIPLIED BY 5*>‘'15 AND THE 

CD Least significant 25 bits are nFfAirM’o and repeal: the 

CO , CONTENTS OF I. THE LEAST SICNIFlCANi 27 RTf- ARE t'!.- 

CO BINED WITH AN EXPONENT AND CONVERTED TO A NoR.'ALIZLD 
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000013 irmirniii 


. 'CD 


FLOATIMG POINT 

NL'NOER, A, milCH HAS A UNIFORM DISTRI3UTI0N 

. CD 
. CD 


ON THE INTERVAL 

from hro to one. 

. CD 

***<!( 

******** 

********** * * * * * Vt * * * * V * 

$(1), 

NBPR. 






ROM* 

LM 

-AO, *0,011 . 

PICK UP A3,5<I). 


OR • 

AO, CON, ,1 . 

SET eiT ZERO ON. 


MI 

A1,C0N1 . 

, GET NEW RANDOM NUH9CR. 


AND 

A2,NASi< . . 

CLEAR BIT 35. 


SA 

A3, *0,011 . 

STORE. 


LA . 

A0,C0N„2 . 

PICK UP CHARACTERISTIC. 


LCF 

AO, IS . 

FLOAT. 


SA 

A1, *1,911 . ■ 

STORE. 


J 

3,B11 . 

RETURN 

CON 

CONI 

MASK 

END 

-1-01 70000001 

+30517578125 

+0377777777777. 



C A 3 II. E!!R0I?S ; I4.0NE 




n 

1 

00 


rmiiji 

■1* 

(10101 

2* 

00101 

3* 

uoitn 

4* 

O'nui 

s* 

01.1101 

6* 

L",’ ! 0 1 

7* 

Oil KM 

8* 


9* 

0UI'’J 

10* 

001U5 

11* 

OO'liM, 

12* 

fjono 

13* 

00 1 1 1 

14* 

00114 

15* 

00116 

16* 

0U121 

17* 

0U12? 

18* 

001 ?i 

19* 

00124 

20* 

00123 

21* 

00 1 26 

22* 

00127 

23* 

001 20 

24* 

UQ131 

25* 


c SUUftOUrifJE RDM IS ENTERED TWICE IN SUCCESION ON 

C ***/<*)•.****** alternate CNlRIbS TO RN2S. RDM USES A FIXED POINT 
C ************ random NUMOER to generate a M0RHAL17.ED FLOATING POINT 

c ************ random numoer having a uniform DISTRIDUTION on the 
c ************ interval from zero to one. a NED VALUE FOR THE riXCD 
C ************ POINT RANDOM NUMDER IS ALSO GENERATED AT EACH ENliiY TO 

C *****tf****** ROH 

DOUBLE PRECISION Q,. SIGMA, X, X1, X2,'X2B 
IITFEGER ORNUM, RNUM, RNUM1 
COMMON/ RANDOM /SIGMA, X2,NSTART,X28,NSTRT0 
1 IF (RNUM.EQ.O) GO TO AO 
RNUm = RNUM 
' IF (X2) 20,60,20 

20 IF(NSTART.EQ.O) GO -TO AO 8USE NON-ZERO X2 AS FIRST R.V. 

IF (RNUMI-ORNUM) 60,40,60 
AO X = X2 
X2 =-0.0 
GO TO GO 

60 CALL RDM(RNUM1,U1) 

CALL R1»M(RNUM1,U2) 

NSTART=--7 
ORNUfl = RNUM1 
•RNUM = RNUM1 

Q = 6.283185307179586A8D0 * U2 


RN7S0001 

RN2S00UJ 

RN2S0UD3 

RN2S000A 

RI-J2S0005 

RNPS0006 

RN2S0007 

RU2SOOiir. 

Ri.'E'iOUOV 

RMESUOIO 

RN2S0011 

RN2S0012 

RN2S0013 

RN2S0014 

RN2S0015 

RN2S0016 

RN2S00I7 

RN2SUU18 

PJ 2 SP 01 ? 

RN2S0020 

RN2S0021 

RN2S0022 

RN2S0023 

RN2S002A 

RN2S0025 


OnOfl'K! 

OiH'lKiO 

OL’OiiOO 

000000 

000000 

OJU-'KiG 

niit.i'it'.n 

Oi-.JOO 

oocyoo 
conooo 
oooooo 
000001 
000003 
000005 
GU0007 
000012 
0000 1 3 
0000 1 5 
000 UI 7 

000 n .?2 

000026 

000030 

000032 

000035 



RIJ2S 

ill, 

DATE 

120178 

OEM 52 

26* 

X1 n DCOS(Q) 

RN2S0026 

000936 

(. 01:3 

2/* 

X2 = OSIN(Q) 

RN2S0027 

000042 

001 54 

28* 

Q = OSORT(-2,DO*DLOG(DBLE(U1))) 

RN2SD028 

0000.46 

00135 

29* ■ 

XI = Q i* X1 

RN2S0029 

000062 

4 

nr.l 

30* 

X2 = (5 * X2 

RN2S0O30 

0015/ 

31 .-1 

X = XI 

RN2S0031 

n0u067 

001 4 U 

3 •* 

80 IF (0ABS(X),.GT. SIGMA) GO TO 1 

RN2S0032 

000071 

00142 

33* 

RETURN 

RN2S0033 

000074 

00143 

34 * 

END 

RN2S003A 

000115 


END OF COMPILATION: NO Ol/'GNOSriCS. 


OF POOR QUAUsY 



